lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -N OuterTrackFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ OuterTrackFinder.java 18 Jul 2007 19:01:06 -0000 1.1
@@ -0,0 +1,1174 @@
+/*
+ * AxialBarrelTrackFinder1.java
+ *
+ * Created on August 18, 2005, 4:17 PM
+ *
+ */
+
+package org.lcsim.contrib.CosminDeaconu;
+
+
+import com.sun.crypto.provider.ARCFOURCipher;
+import org.lcsim.contrib.tracking.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.base.BaseTrackMC;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseTrackerHitMC;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.util.Driver;
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.circle.CircleFitter;
+
+import org.lcsim.contrib.CosminDeaconu.HelicalTrackFit;
+import org.lcsim.contrib.CosminDeaconu.HelicalTrackFitter;
+
+import org.lcsim.geometry.Detector;
+import org.lcsim.detector.converter.compact.DeDetector;
+import org.lcsim.geometry.subdetector.MultiLayerTracker;
+import org.lcsim.geometry.subdetector.DiskTracker;
+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;
+import org.lcsim.contrib.JanStrube.tracking.Helix;
+
+/**
+ * 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 OuterTrackFinder extends Driver
+{
+
+ // Data members
+
+
+ //String arrays are used for easy input here, they are converted to arraylists in the constructor
+ private String[] _input_hit_collections_arr = {"TkrForwardHits","TkrEndcapHits"};
+ //private String[] _input_hit_collections_arr = {"TkrForwardHits","TkrEndcapHits","TkrBarrHits"};
+ private ArrayList<String> _input_hit_collections = null;
+
+
+ //String arrays are used for easy input here, they are converted to arraylists in the constructor
+ private String[] _subdetectors_arr = {"TrackerForward","TrackerEndcap"};
+ //private String[] _subdetectors_arr = {"TrackerForward","TrackerEndcap","TrackerBarrel"};
+ private ArrayList<String> _subdetectors = null;
+
+ 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 = true;
+ private int nAssocMin = 4;
+ private double _min_seed_pt = 0.5; //*Tyler
+ private double _max_phi_sep = (Math.PI)/2; //*Tyler
+
+ private List<TrackerHit> _hits = null;
+ private Vector<List<TrackerHit>> _hitsbylayer = new Vector<List<TrackerHit>>();
+ private List<TrackerHit> _used_hits = new ArrayList<TrackerHit>();
+ private HashMap<TrackerHit,Double> _hit_separation = new HashMap<TrackerHit,Double>();
+
+ private TrackerHitCheater _cheat = new TrackerHitCheater();
+
+
+ private AIDA _aida = AIDA.defaultInstance();
+
+ int _min_layers=4;
+
+ //~Cosmin
+ private HelicalTrackFitter _hfitter = new HelicalTrackFitter();
+ private HelicalTrackFit _hfit = null;
+ private double dz = 100; //
+ private double drphi = 0.0316; //I think this is equivalent to having 1000 for weights on circle fit
+
+ private double[] _h_chisq_dof = {10.,10.};
+
+ private boolean _debug_sysouts = false;
+ private double _debug_sysout_level=6;
+
+ private String _id=""; //this id will(should) be appended to all histograms/events put out.
+ //the purpose is so that the output of multiple AxialBarrelTrackFinder1_1 subdrivers can be distinguished
+ //On a lot of the histograms, this is probably overkill since they'll likely be the same between different subdrivers...
+
+ private boolean _multiple_collections = false; // Will be set to true automatically if subdetectors.size()>1
+ private boolean _useAllLayers = false; //true for barrel, false for forward
+
+ //arrays are used for easy entry, they are converted to arraylists in the constructor
+
+ private Integer[] _usedLayers_arr={0,2,4,6,8,10,12};
+ //private Integer[] _usedLayers_arr={0,2,4,6,8,10,12,14,15,16,17,18};
+ private ArrayList<Integer> _usedLayers = null;
+
+ //arrays are used for easy entry, they are converted to arraylists in the constructor
+ private int [][] _combinedLayers_arr={{4,6}};
+ private ArrayList<ArrayList<Integer>> _combinedLayers = null;
+
+ private boolean _phi_cut = false; //whether or not to use the phicut
+
+ /**
+ * Creates a new instance of AxialBarrelTrackFinder1
+ */
+ public OuterTrackFinder()
+ {
+ //convert arrays to arraylists (so that they can be changed by public methods).
+ _subdetectors = new ArrayList<String>(Arrays.asList(_subdetectors_arr));
+ _input_hit_collections = new ArrayList<String>(Arrays.asList(_input_hit_collections_arr));
+ _usedLayers = new ArrayList<Integer>(Arrays.asList(_usedLayers_arr));
+ _combinedLayers = convert2DArray(_combinedLayers_arr);
+
+ }
+
+ /**HelicalTrackFit
+ * 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];
+
+
+ if (_subdetectors.size()>1) _multiple_collections=true;
+
+ // Get the tracker attributes
+
+ int nlayers=0;
+// double[] radius; //not needed?
+// double[] halflength; //not needed?
+
+ ArrayList subdetector_lengths = new ArrayList(); // this stores the length of each subdetector used. Only makes sense when there are multiple collections
+ ArrayList layer_converter = new ArrayList(); // this arraylist will be used to skip/combine layers, if necessary
+
+ //get nlayers ~Cosmin
+ //--------------------------------------------------------------////
+
+
+ for (String s : _subdetectors)
+ {
+ if (detector.getSubdetector(s).isEndcap()) //for tracker forward or endcap
+ {
+ DiskTracker subdetector = (DiskTracker)detector.getSubdetector(s);
+ subdetector_lengths.add(subdetector.getInnerR().length);
+ nlayers = nlayers + subdetector.getInnerR().length;
+ }
+
+ else if (detector.getSubdetector(s).isBarrel()) // for barrel
+ {
+ MultiLayerTracker subdetector = (MultiLayerTracker)detector.getSubdetector(s);
+ subdetector_lengths.add(subdetector.getInnerR().length);
+ nlayers = nlayers + subdetector.getInnerR().length;
+ }
+ }
+
+
+ //if skipping or combining layers, figure out what to do with each layer ~Cosmin
+ //------------------------------------------------------------------------------//
+ if (!_useAllLayers)
+ {
+ for (int i = 0; i<nlayers; i++) //by default, map each layer to itself
+ {
+ layer_converter.add(i);
+ }
+
+ Set used = new TreeSet();
+
+ for (Integer j : _usedLayers)
+ {
+ int i = j.intValue();
+ used.add(i);
+ if (i>nlayers)
+ {
+ System.out.println("ERROR: THE USED LAYERS ARRAY CONTAINS NON-EXISTING LAYERS");
+ throw (new NullPointerException()); //throw exception if _usedLayers specifies layers that don't exist
+ }
+ }
+
+ //skip layers not in _usedLayers if not using all layers and then map other layers to them
+ Set available_spots = new TreeSet();
+
+ for (int i =0;i<layer_converter.size();i++)
+ {
+ if (!used.contains(i))
+ {
+ layer_converter.set(i,-1); //set unused layers to -1 so we know to skip them
+ available_spots.add(i);
+ }
+ }
+
+ if (_combinedLayers.size()>0) //free up locations held by layers to be combined
+ {
+ for (ArrayList i : _combinedLayers)
+ {
+ for (int j = 1; j<i.size();j++)
+ {
+ used.remove(i.get(j));
+ available_spots.add(i.get(j));
+ }
+ }
+ }
+
+ for (int i = 0; i<layer_converter.size();i++) //map used layers to lowest possible index
+ {
+ if (used.contains(i) && available_spots.size()>0)
+ {
+ for (int j = 0; j<i;j++)
+ {
+ if (available_spots.contains(j))
+ {
+ layer_converter.set(i,j);
+ available_spots.remove(j);
+ available_spots.add(i);
+ break;
+ }
+ }
+ }
+ }
+
+ if (_combinedLayers.size()>0) //combine any layers to be combined
+ {
+ for (ArrayList i : _combinedLayers)
+ {
+ for (int j = 1; j<i.size();j++)
+ {
+ layer_converter.set(((Integer)i.get(j)).intValue(),layer_converter.get(((Integer)(i.get(0))).intValue())); //yuck... why must ArrayLists not work with primitives?
+
+ }
+ }
+ }
+
+ //System.out.println(layer_converter.toString());
+ nlayers = used.size(); //set nlayers to the number of layers actually used
+ }
+
+ gotHere(0);
+
+
+ // Get the SimTrackerHits and metadata
+ //-------------------------------------------------------------------
+
+ _hits = new ArrayList<TrackerHit>();
+ for (String i : _input_hit_collections)
+ {
+ _hits.addAll((List<TrackerHit>)_cheat.makeTrackerHits(event.getSimTrackerHits(i)));
+ }
+
+ // Print out number of hits
+ if (_make_histograms) _aida.cloud1D("nHitsTotal"+_id).fill(_hits.size());
+
+
+ // Create/clear vector of TrackerHits for each layer
+
+ if (_hitsbylayer.size()==0)
+ {
+ for (int layer = 0; layer < nlayers; layer++)
+ {
+ _hitsbylayer.add(layer, new Vector<TrackerHit>());
+ }
+ }
+ else
+ {
+ for (int layer = 0; layer < nlayers; layer++)
+ {
+ _hitsbylayer.get(layer).clear();
+ }
+ }
+
+ for (TrackerHit hit : _hits)
+ {
+
+ int layer = ((BaseTrackerHitMC)hit).getSimHits().get(0).getLayer();
+
+ if (_multiple_collections) //Stack layers if multiple hit collections so we don't have a bunch of layer 0's and such'
+ {
+ for (int i = 1; i<_subdetectors.size();i++)
+ {
+ if (((BaseTrackerHitMC)hit).getSimHits().get(0).getSubdetector().getName().equals(_subdetectors.get(i)))
+ {
+ int offset = 0;
+ for (int j = i-1; j>-1;j--)
+ {
+ int this_offset = ((Integer)subdetector_lengths.get(j)).intValue();
+ offset=offset+this_offset;
+ }
+// System.out.println("offset:"+offset);
+ layer = layer+offset;
+ }
+ }
+ }
+
+ if (!_useAllLayers) //if not using all layers, go through the necessary hoopla
+ {
+ if (((Integer)layer_converter.get(layer)).intValue()==-1) //skip skipped layers (which have a sentinel value of -1)
+ {
+ continue;
+ }
+
+ layer = ((Integer)layer_converter.get(layer)).intValue(); // map layers to correct order so that the layer collection only has length of nlayers
+
+ }
+// System.out.println(hit.getSubdetector().getName()+"layer: "+layer); //code for debugging
+
+ _hitsbylayer.elementAt(layer).add(hit);
+ }
+
+ gotHere(1);
+
+ // Create map from hits to separation from nearest hit in layer
+ _hit_separation.clear();
+ for (int layer = 0; layer < nlayers; layer++)
+ {
+ List<TrackerHit> hits = _hitsbylayer.get(layer);
+ for (TrackerHit 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>();
+ HashSet<MCParticle> enough_layer_particles = new HashSet<MCParticle>();
+
+ int ntracks_all_layers = 0;
+ int ntracks_enough_layers=0;
+ for (MCParticle mc_particle : mc_particles)
+ {
+ boolean layers_hit[] = new boolean[nlayers];
+ for (TrackerHit hit : _hits)
+ {
+
+ int layer = ((BaseTrackerHitMC)hit).getSimHits().get(0).getLayer();
+
+
+ if (_multiple_collections) //Stack layers if multiple hit collections so we don't have a bunch of layer 0's and such'
+ {
+ for (int i = 1; i<_subdetectors.size();i++)
+ {
+ if (((BaseTrackerHitMC)hit).getSimHits().get(0).getSubdetector().getName().equals(_subdetectors.get(i)))
+ {
+ int offset = 0;
+ for (int j = i-1; j>-1;j--)
+ {
+ int this_offset = ((Integer)subdetector_lengths.get(j)).intValue();
+ offset=offset+this_offset;
+ }
+ // System.out.println("offset:"+offset);
+ layer = layer+offset;
+ }
+ }
+ }
+
+ if (!_useAllLayers) //if not using all layers, go through the necessary hoopla
+ {
+ if (((Integer)layer_converter.get(layer)).intValue()==-1) //skip skipped layers (which have a sentinel value of -1)
+ {
+ continue;
+ }
+
+ layer = ((Integer)layer_converter.get(layer)).intValue(); // map layers to correct order so that the layer collection only has length of nlayers
+
+ }
+ // System.out.println(hit.getSubdetector().getName()+"layer: "+layer); //code for debugging
+
+ if(_make_histograms) _aida.cloud1D("Number of MC Particles per hit").fill(((BaseTrackerHitMC)hit).mcParticles().size());
+
+ if (mc_particle == ((BaseTrackerHitMC)hit).getSimHits().get(0).getMCParticle())
+ {
+ layers_hit[layer] = true;
+ }
+ }
+
+
+
+ int num_layers = 0;
+ for (int layer =0; layer<nlayers;layer++)
+ {
+ if (layers_hit[layer])
+ {
+ num_layers++;
+ }
+ }
+
+ if (num_layers >= _min_layers)
+ {
+ ntracks_enough_layers++;
+ enough_layer_particles.add(mc_particle);
+ if (_make_histograms) _aida.cloud1D(">"+_min_layers+"-hit MCParticle momentum"+_id).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;
+
+ gotHere(2);
+
+ // Create list of tracks
+ //--------------------------------------------------------------------------------------------//
+ //ArrayList<StandaloneAxialBarrelTrack1> tracklist = new ArrayList<StandaloneAxialBarrelTrack1>();
+ ArrayList<StandaloneOuterTrack> tracklist = new ArrayList<StandaloneOuterTrack>();
+ ArrayList<StandaloneAxialBarrelTrack1> faketracks = new ArrayList<StandaloneAxialBarrelTrack1>(); //IS THIS EVER USED?
+ // 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>();
+ gotHere(4.900001);
+ List<TrackerHit> hits1 = _hitsbylayer.get(layer1);
+ List<TrackerHit> hits2 = _hitsbylayer.get(layer2);
+ List<TrackerHit> hits3 = _hitsbylayer.get(layer3);
+
+ double[] x = new double[nlayers];
+ double[] y = new double[nlayers];
+ double[] z = new double[nlayers];
+ //double[] weight = new double[nlayers];
+
+ // Try all combinations of hits
+ for (TrackerHit 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.getPosition()[2] > 0.0); // find +/- z
+ //double zrel = (hit1.getPoint()[2]/radius[layer1]);
+ x[0] = hit1.getPosition()[0];
+ y[0] = hit1.getPosition()[1];
+ z[0] = hit1.getPosition()[2];
+ // weight[0] = 1000.0; // set to achieve reasonable chisq
+ for (TrackerHit 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 ((hit2.getPosition()[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.getPosition()[0];
+ y[1] = hit2.getPosition()[1];
+ z[1] = hit2.getPosition()[2];
+ //weight[1] = 1000.0; // set to achieve reasonable chisq
+ for (TrackerHit hit3 : hits3)
+ {
+ gotHere(4.950001);
+ 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 ((hit3.getPosition()[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.getPosition()[0];
+ y[2] = hit3.getPosition()[1];
+ z[2] = hit3.getPosition()[2];
+ // weight[2] = 1000.0; // set to achieve reasonable chisq
+
+ // Perform circle fit with these three points, marking layers as used
+ int nhits = 3;
+ Vector<TrackerHit> hitset = new Vector<TrackerHit>();
+ 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);
+ boolean fitok = _hfitter.fit(x,y,z,get_drphi(nhits),get_dz(nhits),nhits);
+
+ if (!fitok) System.out.println("Warning: seed fit failed!");
+ //_fit = _fitter.getfit();
+ _hfit = _hfitter.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());
+ double pc = (0.00029979*b_field)/(_hfit.parameters()[2]);
+ //if (Math.abs(_fit.dca()) <_seed_ip_dca && Math.abs(pc) > _min_seed_pt )
+ gotHere(4.990001);
+ _aida.cloud1D("LineChisQ").fill(_hfit.chisq()[1]);
+ _aida.cloud1D("CircChisQ").fill(_hfit.chisq()[0]);
+ if (Math.abs(_hfit.parameters()[0]) < _seed_ip_dca && Math.abs(pc) > _min_seed_pt && _hfit.chisq()[1]/3.0<_h_chisq_dof[1])
+ {
+ //System.out.println("Passed chisq cut");
+ //_aida.cloud1D("CircleChisQ").fill(_hfit.chisq()[0]);
+ gotHere(5);
+ // Fill seed histograms
+ if (_make_histograms)
+ {
+ //_aida.cloud1D("dca_seed"+_id).fill(_fit.dca());
+ _aida.cloud1D("dca_seed"+_id).fill(_hfit.parameters()[0]);
+
+ //_aida.cloud1D("curvature_seed"+_id).fill(_fit.curvature());
+ _aida.cloud1D("curvature_seed"+_id).fill(_hfit.parameters()[2]);
+ }
+
+ for (int layer = nlayers-1; layer >= 0; layer--)
+ {
+ gotHere(5.1);
+ if (used_layers.contains(layer)) continue;
+ gotHere(5.2);
+ // Add nearest hit if unambiguous -- check order of logic here
+ List<TrackerHit> hits = _hitsbylayer.get(layer);
+ //System.out.println("hits size: "+hits.size());
+ double nearest_dist = 1000000.0;
+ double next_nearest_dist = 1000000.0;
+ TrackerHit nearest_hit = null;
+
+ for (TrackerHit hit : hits)
+ {
+ gotHere(5.29);
+ if (_used_hits.contains(hit)) continue; // skip used hits
+ if ((hit.getPosition()[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)
+ gotHere(5.295);
+ Helix h = _hfit.getHelix();
+ BasicHep3Vector p = new BasicHep3Vector(hit.getPosition());
+ gotHere(5.296);
+ double distance = Math.abs(h.getSignedClosestDifferenceToPoint(new SpacePoint(p)));
+ //System.out.println("distance: "+distance);
+ gotHere(5.3);
+ if (distance < nearest_dist)
+ {
+ gotHere(5.4);
+ next_nearest_dist = nearest_dist;
+ //nearest_dist = Math.abs(_fit.dca());
+ nearest_dist = distance;
+
+ nearest_hit = hit;
+ }
+ }
+
+ // Make histograms of nearest hits
+ if (_make_histograms)
+ {
+ if (nearest_dist < 10.0)
+ {
+ _aida.cloud1D("dca_add1_nearest"+layer).fill(nearest_dist);
+ }
+ if (next_nearest_dist < 10.0)
+ {
+ _aida.cloud1D("dca_add1_nextnearest"+layer).fill(next_nearest_dist);
+ }
+ }
+
+ // require nearby hit that is well separated
+ if (nearest_dist < _pass1_hit_dca && next_nearest_dist > _pass1_hit_isolation)
+ {
+ gotHere(6);
+ used_layers.add(layer); // Mark layer used
+ hitset.setElementAt(nearest_hit,layer);
+ x[nhits] = nearest_hit.getPosition()[0];
+ y[nhits] = nearest_hit.getPosition()[1];
+ z[nhits] = nearest_hit.getPosition()[2];
+ //weight[nhits] = 1000.0; // set to achieve reasonable chisq
+ nhits++;
+
+ // Refit after adding each hit
+ //fitok = _fitter.fit(x,y, weight, nhits);
+ fitok = _hfitter.fit(x,y,z,get_drphi(nhits),get_dz(nhits),nhits);
+ if (!fitok) System.out.println("Warning: intermediate fit failed!");
+ }
+ }
+
+ gotHere(6.5);
+
+ // 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--)
+ {
+ gotHere(6.6);
+ if (used_layers.contains(layer)) continue;
+
+ // Find nearest hit
+ List<TrackerHit> hits = _hitsbylayer.get(layer);
+ double nearest_dist = 1000000.0; // 1km
+ double next_nearest_dist = 1000000.0;
+ TrackerHit nearest_hit = null;
+ gotHere(6.8);
+ for (TrackerHit hit : hits)
+ {
+ if (_used_hits.contains(hit)) continue; // skip used hits
+ if ((hit.getPosition()[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)
+ gotHere(6.9);
+ double distance = Math.abs(_hfit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(hit.getPosition()))));
+ if( distance < nearest_dist)
+ {
+ gotHere(7);
+ next_nearest_dist = nearest_dist;
+ //nearest_dist = Math.abs(_fit.dca())
+ nearest_dist = distance;
+ nearest_hit = hit;
+ }
+ }
+ // Make histograms of nearest hits
+ if (_make_histograms)
+ {
+ if (nearest_dist < 10.0)
+ {
+ _aida.cloud1D("dca_add2_nearest"+_id).fill(nearest_dist);
+ }
+ if (next_nearest_dist < 10.0)
+ {
+ _aida.cloud1D("dca_add2_nextnearest"+_id).fill(next_nearest_dist);
+ }
+ }
+
+ gotHere(7.5);
+ // require nearby hit
+ if (nearest_dist < _pass2_hit_dca)
+ {
+ gotHere(8);
+ used_layers.add(layer); // Mark layer used
+ hitset.setElementAt(nearest_hit,layer);
+ x[nhits] = nearest_hit.getPosition()[0];
+ y[nhits] = nearest_hit.getPosition()[1];
+ z[nhits] = nearest_hit.getPosition()[2];
+ //weight[nhits] = 1000.0; // set to achieve reasonable chisq
+ nhits++;
+ }
+ }
+ } // if less than 5 layers
+
+ // Reset reference position
+ //_fitter.setreferenceposition(0.0,0.0);
+ gotHere(9);
+ // If have successfully added hits, perform final fit
+ if (used_layers.size() >= _min_layers)
+ {
+ //fitok = _fitter.fit(x,y, weight, nhits);
+ fitok = _hfitter.fit(x,y,z,get_drphi(nhits),get_dz(nhits),nhits);
+ if (!fitok) System.out.println("Warning: final fit failed!");
+
+ //_fit = _fitter.getfit();
+ _hfit = _hfitter.getFit();
+ double[] phivalues2 = new double[nlayers];//~Cosmin, changed a 5->nlayers
+ int k2 = 0;
+
+ //Create Array of Phi-hit values *Tyler
+ for (TrackerHit hit : hitset)
+ {
+ if(hit == null) continue;
+ k2++;
+ double[] hitpoint2 = hit.getPosition();
+ 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[nlayers-1] == 0)
+ {
+ phivalues2[nlayers-1] = phivalues2[0];
+ }
+
+ int b2=0;
+ for(int ps2=0; ps2<nlayers; ps2++)
+ {
+
+ //Test hits to be sure all hits have a reasonable difference in phi *Tyler
+ for(int ts2=0; ts2<nlayers; ts2++)
+ {
+ double phidifference2 = Math.abs(phivalues2[ps2]-phivalues2[ts2]);
+ _aida.cloud1D("PhiDifference").fill(phidifference2);
+ 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 || !_phi_cut)
+ {
+
+ // If fit is good enough, fill final histograms
+ //if (_fit.chisq()/nhits < _chisq_dof)
+ _aida.cloud1D("ChisqL before check").fill(_hfit.chisq()[1]);
+ _aida.cloud1D("ChisqC before check").fill(_hfit.chisq()[0]);
+ if (_hfit.chisq()[0]/nhits < _h_chisq_dof[0] && _hfit.chisq()[1]/nhits < _h_chisq_dof[1])
+ {
+ // Mark hits as used - find all MCParticles that contribute
+ HashSet<MCParticle> mc_particleset = new HashSet<MCParticle>();
+ int n = 0;
+ for (TrackerHit 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(((BaseTrackerHitMC)hit).getSimHits().get(0).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 (TrackerHit hit : hitset)
+ {
+ if (hit == null) continue;
+ if (particle == ((BaseTrackerHitMC)hit).getSimHits().get(0).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)
+ if (enough_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"+_id).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"+_id).fill(_fit.chisq());
+ _aida.cloud1D("chisq[0]_pass"+_id).fill(_hfit.chisq()[0]);
+ _aida.cloud1D("chisq[1]_pass"+_id).fill(_hfit.chisq()[1]);
+ //_aida.cloud1D("dca_pass"+_id).fill(_fit.dca());
+ _aida.cloud1D("dca_pass"+_id).fill(_hfit.parameters()[0]);
+ _aida.cloud1D("nhits_pass"+_id).fill(used_layers.size());
+ _aida.cloud1D("purity"+_id).fill((double)nhits_max/(double)nhits);
+ }
+
+ // Create StandaloneAxialBarrelTrack1
+ //StandaloneAxialBarrelTrack1 track = new StandaloneAxialBarrelTrack1(b_field, _module_length, _fit, hitset, majority_particle, purity);
+ StandaloneOuterTrack track = new StandaloneOuterTrack(b_field, _module_length, _hfit, 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"+_id).fill(nfound);
+ //if (ntracks_all_layers > 0)
+ if (ntracks_enough_layers > 0)
+ {
+ //_aida.cloud1D("efficiency"+_id).fill((double)nfound/(double)ntracks_all_layers);
+ _aida.cloud1D("efficiency"+_id).fill((double)nfound/(double)ntracks_enough_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<TrackerHit> unused_hits = _hits;
+ int i = 0;
+ for (TrackerHit hit : _used_hits)
+ {
+ i++;
+ unused_hits.remove(hit);
+ }
+
+ event.put(_output_hit_collection, unused_hits, TrackerHit.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"+_id).isConverted())
+ {
+ _aida.cloud1D("Found MCParticle momentum"+_id).convert(50, 0.0, 50.0);
+ }
+ if (!_aida.cloud1D("5-hit MCParticle momentum"+_id).isConverted())
+ {
+ _aida.cloud1D("5-hit MCParticle momentum"+_id).convert(50, 0.0, 50.0);
+ }
+
+ hfactory.divide("efficiency vs. momentum"+_id,
+ _aida.cloud1D("Found MCParticle momentum"+_id).histogram(),
+ _aida.cloud1D("5-hit MCParticle momentum"+_id).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 = "XXXX")
+ *
+ * @param input_hit_collection Name of input SimTrackerHit collection
+ *
+ */
+ public void setInputHits(String input_hit_collection)
+ {
+ _input_hit_collections.clear();
+ _input_hit_collections.add(input_hit_collection);
+ }
+
+ /**
+ * Set name of SimTrackerHit collection to make tracks from (default = "XXXXX")
+ *
+ * @param input_hit_collection String array of names of input SimTrackerHit collection
+ *
+ */
+ public void setInputHits(String[] input_hit_collection)
+ {
+ _input_hit_collections.clear();
+ _input_hit_collections = new ArrayList<String>(Arrays.asList(input_hit_collection));
+ }
+
+
+ /**
+ * Set name of subdetectors to get layer info from (default = "XXXXX")
+ *
+ * @param subdetectors String array of names of subdetectors
+ *
+ */
+ public void setSubdetectors(String[] subdetectors)
+ {
+ _subdetectors.clear();
+ _subdetectors = new ArrayList<String>(Arrays.asList(subdetectors));
+ }
+
+ /**
+ * Sets whether or not all layers will be used in a subdetector (default = "XXXXX")
+ *
+ * @param useAllLayers true to use all, false otherwise
+ *
+ */
+ public void setUseAllLayers(boolean useAllLayers)
+ {_useAllLayers = useAllLayers;}
+
+ /**
+ * Sets the layers to be used by the finder (default = "XXXXX")
+ *
+ * @param usedLayers int array of layers to be used
+ *
+ */
+ public void setUsedLayers(int[] usedLayers)
+ {
+ _usedLayers.clear();
+ _usedLayers = new ArrayList(Arrays.asList((usedLayers)));
+ }
+
+ public void setCombinedLayers(int[][] combinedLayers)
+ {
+ _combinedLayers.clear();
+ _combinedLayers = convert2DArray(combinedLayers);
+ }
+
+ /**
+ * Set name of SimTrackerHit collection for unused hits (default = "" == not written out)
+ *
+ * @param output_hit_collection Name of output collection for unused SimTrackerHits
+ *if (ntracks_all_layers > 0)
+ */
+ 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)
+ *
+ */
[truncated at 1000 lines; 178 more skipped]
lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -N HelicalTrackFitter.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelicalTrackFitter.java 18 Jul 2007 19:01:06 -0000 1.1
@@ -0,0 +1,147 @@
+/*
+ * HelicalTrackFitter.java
+ *
+ * Created on March 25, 2006, 6:11 PM
+ *
+ * $Id: HelicalTrackFitter.java,v 1.1 2007/07/18 19:01:06 cozzyd Exp $
+ */
+
+package org.lcsim.contrib.CosminDeaconu;
+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.line.SlopeInterceptLineFit;
+import org.lcsim.fit.line.SlopeInterceptLineFitter;
+/**
+ *
+ * @author Norman Graf
+ */
+public class HelicalTrackFitter
+{
+ private CircleFitter _cfitter = new CircleFitter();
+ private SlopeInterceptLineFitter _lfitter = new SlopeInterceptLineFitter();
+
+ private HelicalTrackFit _fit;
+
+ private double[] _circleParameters = new double[3];
+
+ /** Creates a new instance of HelicalTrackFitter */
+ public HelicalTrackFitter()
+ {
+ }
+
+ public boolean fit(double[] x, double[] y, double[] z, double[] drphi, double[] dz, int np)
+ {
+ // fit a circle in the x-y plane
+ double[] wrphi = new double[np];
+ for(int i=0; i<np; ++i)
+ {
+ wrphi[i] = 1/(drphi[i]*drphi[i]);
+ }
+ boolean success = _cfitter.fit(x, y, wrphi, np);
+ 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
+ //TODO check whether hit has z information
+ // TODO double-check path length calculation
+ double[] s = new double[np];
+ double[] wz = new double[np];
+ for(int i=0; i<np; ++i)
+ {
+ s[i] = s(radius, xc, yc, x[i], y[i], 0., 0.);
+ wz[i] = 1/(dz[i]*dz[i]);
+ }
+
+ success = _lfitter.fit(s, z, dz, np);
+ if(!success) return false;
+ SlopeInterceptLineFit lfit = _lfitter.getFit();
+
+ // TODO convert fit results to track parameters and covariance matrix.
+ // now have the fits, need to assemble the track parameters...
+
+ // see org.lcsim.event.Track.Parameters
+ // The track parameters for LCD are defined as follows
+ // <table>
+ // <tr><th>Index</th><th>Meaning</th></tr>
+ // <tr><td> 0 </td><td> d0 = XY impact parameter </td><tr>
+ // <tr><td> 1 </td><td> phi0 </td><tr> </td><tr>
+ // <tr><td> 2 </td><td> omega = 1/curv.radius (negative for negative tracks) </td><tr>
+ // <tr><td> 3 </td><td> z0 = z of track origin (z impact parameter) </td><tr>
+ // <tr><td> 4 </td><td> s = tan lambda </td><tr>
+ // </table>
+
+ // tanlambda = line fit slope
+ // z0 = line fit intercept
+ //TODO should shift s-z line fit to centroid of x-y fit and introduce covariance terms to be correct.
+
+
+// System.out.println("circle fit: "+cfit);
+// System.out.println("line fit: "+lfit);;
+ 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();
+ //System.out.println("chisqL:"+chisq[1]);
+ ndf[1] = lfit.ndf();
+ // TODO fix this so CircleFit returns ndf
+ ndf[0] = lfit.ndf();
+
+ // d0, xy impact parameter.
+ 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();
+
+ // TODO fill in covariance matrix
+// for(int i=0; i<5; ++i)
+// {
+// System.out.println("pars["+i+"]= "+pars[i]);
+// }
+ _fit = new HelicalTrackFit(pars, cov, chisq, ndf);
+ _fit.setCircleParameters(_circleParameters);
+ return true;
+ }
+
+ public HelicalTrackFit getFit()
+ {
+ return _fit;
+ }
+
+ 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);
+ }
+
+}
\ No newline at end of file
lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -N StandaloneOuterTrack.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ StandaloneOuterTrack.java 18 Jul 2007 19:01:06 -0000 1.1
@@ -0,0 +1,524 @@
+/*
+ * StandaloneAxialBarrelTrack1.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.
+ */
+
+package org.lcsim.contrib.CosminDeaconu;
+
+import org.lcsim.contrib.tracking.*;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.event.MCParticle;
+import java.util.*;
+import org.lcsim.contrib.CosminDeaconu.HelicalTrackFit;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+
+/**
+ * Implements the Track interface for tracks found using only outer tracker information by TrackFinder
+ *
+ * @author tknelson. Forked by Cosmin Deaconu.
+ * @see AxialBarrelTrackFinder
+ *
+ */
+public class StandaloneOuterTrack implements Track
+{
+
+ // Data members
+ boolean _fit_success = false;
+ double _chi2 = 0.0;
+ double[] _h_chi2 = {0.0,0.0};
+ int _ndof = 0;
+ double _dedx = 0.0;
+ double _dedx_error = 0.0;
+ SymmetricMatrix _error_matrix = new SymmetricMatrix(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;
+ MCParticle _majority_particle;
+ double _purity = 0.0;
+
+ private boolean _helical_fit=false;
+ // Static
+ static int _type = 2; // ???
+
+
+ // Constructors
+ //=============
+
+ /**
+ * Create a default instance of StandaloneAxialBarrelTrack1
+ */
+ public StandaloneOuterTrack()
+ {
+ }
+
+ /**
+ * 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
+ */
+
+
+ //~Cosmin
+ public StandaloneOuterTrack(double b_field, double module_length, HelicalTrackFit fit, List<TrackerHit> hits, MCParticle majority_particle, double purity)
+ {
+ _b_field = b_field;
+ _module_length = module_length;
+
+ _fit_success = true;
+ _h_chi2 = fit.chisq();
+ _helical_fit=true;
+ _ndof = hits.size();
+
+ _reference_point[0] = 0.; //yeah, this is bad
+ _reference_point[1] = 0.; //yeah, this is bad ~/Cosmin
+ _reference_point[2] = 0.0;
+
+ _tracker_hits = hits;
+ double[] _ip = {0.,0.,0.};
+ _track_parameters[0] = fit.parameters()[0];
+ _track_parameters[1] = fit.parameters()[1]-Math.PI; // who ordered this?
+ _track_parameters[2] = fit.parameters()[2]; // convert 1/mm to 1/cm
+ _track_parameters[3] = fit.parameters()[3];
+ _track_parameters[4] = fit.parameters()[4];
+
+ _isreferencepoint_pca = true;
+
+ _majority_particle = majority_particle;
+ _purity = purity;
+
+ _tracks.add(this);
+
+ }
+
+
+
+ // Methods needed to further initialize a StandaloneAxialBarrelTrack
+ //=======================================================
+
+ /**
+ * 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 error matrix of the track fit
+ * <p>
+ * <b>N.B. Currently returns an empty matrix</b>
+ */
+ public SymmetricMatrix getErrorMatrix()
+ {
+ return _error_matrix;
+ }
+
+ /**
+ * Retreive the chisquared of the track fit
+ */
+ public double getChi2()
+ {
+ if(_helical_fit)
+ {
+ _chi2=(_h_chi2[0]+_h_chi2[1])/2.0;
+
+ }
+ 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)
+ {
+ if(hit != null)
+ {
+ 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;
+ }
+ /**
+ *Get majority MCParticle associated with track
+ */
+ public MCParticle getMCParticle()
+ {
+ return _majority_particle;
+ }
+
+ /**
+ *Return percentage of hits from majority particle
+ */
+ public double purity()
+ {
+ return _purity;
+ }
+
+ public double chiSquared()
+ {
+ double resolution = .007; //arclength resolution in mm
+
+ //Get track parameters
+ double dca = _track_parameters[0];
+ double phi = _track_parameters[1];
+
+ // double rtrack = -1.0/_track_parameters[2]; //radius of track in mm, I hope // **GOLD**
+
+ double rtrack = -1.0/_track_parameters[2];
+
+ // int sign = (int)Math.signum(rtrack); // **GOLD**
+
+ int sign = (int)Math.signum(rtrack);
+
+ double charge = getCharge();
+ double tp2 = _track_parameters[2];
+
+ //System.out.println( "Sign is: " +sign+ " and rtrack is: " +rtrack);
+ //System.out.println( "TP2: " +tp2+ "Charge: " +charge);
+
+ /* GOLD
+ //Calculate center of circle defined by track, in polar coordinates
+ double rcenter = Math.abs(rtrack) - (sign*dca);
+ double phicenter = phi-sign*0.5*Math.PI;
+ END GOLD! */
+
+ double rcenter = Math.abs(rtrack) - (sign*dca);
+ double phicenter = sign*phi-0.5*Math.PI;
+
+ /*
+ if (sign > 0){
+ // phicenter = phi - (sign*0.5*Math.PI); // **GOLD**
+
+ //phicenter = sign*phi - (0.5*Math.PI);
+
+ }
+ if (sign < 0){
+ // phicenter = -1*phi+(sign*0.5*Math.PI); //**GOLD**
+
+ //phicenter = sign*phi-(0.5*Math.PI);
+
+ }
+ */
+
+ double wssr = 0.0; //sum of squares of residuals
+
+ //Calculate residuals
+ for(SimTrackerHit hit : _sim_tracker_hits)
+ {
+ if(hit != null)
+ {
+ //Find polar coordinates of hit
+ double rhit = Math.sqrt(hit.getPoint()[0]*hit.getPoint()[0] + hit.getPoint()[1]*hit.getPoint()[1]);
+ double phihit = Math.atan2(hit.getPoint()[1],hit.getPoint()[0]);
+ if(phihit > Math.PI) phihit -= 2*Math.PI;
+ if(phihit < -1.0 * Math.PI) phihit += 2*Math.PI;
+
+ //Find phi coordinate of the point on the track with same radius as hit.
+ //This equation comes from solving for the phi coordinate of the intersection of two circles.
+ //There are, of course, two solutions. I hope I did the right thing to pick which one is good.
+ double deltaPhi = Math.acos((rhit*rhit + rcenter*rcenter - rtrack*rtrack)/(2.0 * rhit * rcenter));
+ double phitrack = 0.0;
+ double phitrack1 = sign*phicenter + deltaPhi;
+ double phitrack2 = sign*phicenter - deltaPhi;
+ if(Math.abs(phihit - phitrack1) < Math.abs(phihit - phitrack2)) phitrack = phitrack1;
+ else phitrack = phitrack2;
+ if(Math.abs(phitrack) > Math.PI)
+ {
+ if(phitrack>0) phitrack -= 2*Math.PI;
+ if(phitrack<0) phitrack += 2*Math.PI;
+ }
+ double dPhi = phihit - phitrack;
+ //Find weighted squared residual
+ double residual = (dPhi * dPhi) * (rhit * rhit)/(resolution * resolution);
+ wssr += residual;
+ //System.out.println("R Hit is: " +rhit+ " dPhi is: " +dPhi);
+ //double logchisquared = Math.log10(wssr);
+ //aida.cloud2D("Log ChiSquared vs. Phi").fill(logchisquared, phi);
+ }
+
+
+ }
+ return wssr;
+ }
+}
\ No newline at end of file
lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -N HelicalTrackFit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelicalTrackFit.java 18 Jul 2007 19:01:07 -0000 1.1
@@ -0,0 +1,101 @@
+/*
+ * HelicalTrackFit.java
+ *
+ * Created on March 25, 2006, 6:11 PM
+ *
+ * $Id: HelicalTrackFit.java,v 1.1 2007/07/18 19:01:07 cozzyd Exp $
+ */
+
+package org.lcsim.contrib.CosminDeaconu;
+
+import com.sun.crypto.provider.ARCFOURCipher;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.contrib.JanStrube.tracking.Helix; //~Cosmin
+
+/**
+ * Represents the result of a fit to a helical track.
+ * @author Norman Graf
+ */
+public class HelicalTrackFit
+{
+ // fit independently to a circle in x-y and a line in s-z
+ // first is circle, second is line
+ private double[] _chisq = new double[2];
+ private int[] _ndf = new int[2];
+ private double[] _parameters;
+ private SymmetricMatrix _covmatrix;
+ private double[] _circleParameters = new double[3];
+
+ /** Creates a new instance of HelicalTrackFit */
+ public HelicalTrackFit(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf)
+ {
+ _parameters = pars;
+ _covmatrix = cov;
+ _chisq = chisq;
+ _ndf = ndf;
+ }
+
+ public double[] parameters()
+ {
+ return _parameters;
+ }
+ public SymmetricMatrix covariance()
+ {
+ return _covmatrix;
+ }
+
+ public double[] chisq()
+ {
+ return _chisq;
+ }
+
+ public int[] ndf()
+ {
+ return _ndf;
+ }
+
+ public String toString()
+ {
+ StringBuffer sb = new StringBuffer("HelicalTrackFit: \n");
+
+ sb.append("d0= "+_parameters[0]+"\n");
+ sb.append("phi0= "+_parameters[1]+"\n");
+ sb.append("curvature: "+_parameters[2]+"\n");
+ sb.append("z0= "+_parameters[3]+"\n");
+ sb.append("tanLambda= "+_parameters[4]+"\n");
+
+ return sb.toString();
+ }
+
+ public void setCircleParameters(double[] pars)
+ {
+ System.arraycopy(pars,0,_circleParameters,0,3);
+ }
+
+ public double[] circleParameters()
+ {
+ return _circleParameters;
+ }
+
+ /** This function returns the Helix associated with the track. The Helix class used is in org.lcsim.contribum.Jan
+ * Someone should really check my math, since I'm probably confused. ~Cosmin
+ */
+ public Helix getHelix()
+ {
+
+ double x=-_parameters[0]*Math.sin(_parameters[1]);
+ double y=_parameters[0]*Math.cos(_parameters[1]);
+ double z= _parameters[3];
+ double phi=_parameters[1];
+ double C = 1.0/_parameters[2];
+ double lambda = Math.atan(_parameters[4]);
+
+
+ return new Helix(new SpacePoint(new BasicHep3Vector(x,y,z)),C,phi,lambda);
+
+ }
+
+}