Commit in hps-java/src/main on MAIN
java/org/lcsim/hps/users/mgraham/MainJASDriver.java+31added 1.1
                                /FastTrackAnalysisDriver.java+178-891.2 -> 1.3
                                /DetailedAnalysisDriver.java+68-641.2 -> 1.3
java/org/lcsim/hps/recon/tracking/FindableTrack.java+1-11.1 -> 1.2
resources/org/lcsim/hps/steering/HPSTrackingBunchedMC.lcsim+58added 1.1
java/org/lcsim/hps/recon/vertexing/BilliorVertex.java+17-141.2 -> 1.3
java/org/lcsim/hps/util/ConvertToStdhep.java+310-1201.4 -> 1.5
resources/org/lcsim/hps/recon/tracking/strategies/HPS-Test-Lyr50.xml+30added 1.1
+693-288
3 added + 5 modified, total 8 files
Some new stuff/changes for running vertexing simulation again.

hps-java/src/main/java/org/lcsim/hps/users/mgraham
MainJASDriver.java added at 1.1
diff -N MainJASDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MainJASDriver.java	29 Aug 2012 15:40:46 -0000	1.1
@@ -0,0 +1,31 @@
+package org.lcsim.hps.users.mgraham;
+
+import org.lcsim.hps.recon.tracking.HPSSVTSensorSetup;
+import org.lcsim.hps.recon.tracking.TrackerDigiDriver;
+import org.lcsim.hps.recon.tracking.TrackerReconDriver;
+import org.lcsim.util.Driver;
+
+/**
+ * Driver for track reconstruction and analysis of HPS detector for execution in JAS.
+ *
+ * @author M. Graham 
+ */
+public final class MainJASDriver extends Driver {
+
+
+
+    public MainJASDriver() {
+        add(new HPSSVTSensorSetup());
+        add(new TrackerDigiDriver());
+//        add(new TrackerReconDriver());
+        TrackerReconDriver trd=new TrackerReconDriver();
+        trd.setStripMaxSeparation(20.0);
+        trd.setStripTolerance(1.0);
+        trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-Test-4pt0.xml");
+        add(trd);
+        add(new DetailedAnalysisDriver(10));
+       
+    }
+
+  
+}

hps-java/src/main/java/org/lcsim/hps/users/mgraham
FastTrackAnalysisDriver.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- FastTrackAnalysisDriver.java	1 Aug 2011 18:25:20 -0000	1.2
+++ FastTrackAnalysisDriver.java	29 Aug 2012 15:40:46 -0000	1.3
@@ -48,33 +48,52 @@
 public class FastTrackAnalysisDriver extends Driver {
 
     String[] detNames = {"Tracker"};
-    Integer _minLayers=8;
-    Integer[] nlayers = {8};
+    Integer _minLayers = 8;
+    Integer[] nlayers = {10};
     int nevt = 0;
     double xref = 50.0; //mm
     public String outputTextName = "myevents.txt";
     FileWriter fw;
     PrintWriter pw;
     double[] beamsize = {0.001, 0.02, 0.02};
-
-    String _config="3pt4";
+    String _config = "Test";
 // flipSign is a kludge...
 //  HelicalTrackFitter doesn't deal with B-fields in -ive Z correctly
 //  so we set the B-field in +iveZ and flip signs of fitted tracks
 //  note:  this should be -1 for Test configurations and +1 for Full (v3.X and lower) configurations
 //  this is set by the _config variable (detType in HeavyPhotonDriver)
-      int flipSign=1;
+    int flipSign = 1;
+
     public FastTrackAnalysisDriver(int trackerLayers, int mintrkLayers, String config) {
         nlayers[0] = trackerLayers;
-        _minLayers=mintrkLayers;
-        _config=config;
-        if(_config.contains("Test")) flipSign=-1;
+        _minLayers = mintrkLayers;
+        _config = config;
+        if (_config.contains("Test")) {
+            flipSign = -1;
+        }
+    }
+
+    public FastTrackAnalysisDriver() {
+    }
+
+    public void setNLayers(int trackerLayers) {
+        nlayers[0] = trackerLayers;
+    }
+
+    public void setMinLayers(int mintrkLayers) {
+        _minLayers = mintrkLayers;
+    }
+
+    public void setConfig(String config) {
+        _config = config;
+        if (_config.contains("Test")) {
+            flipSign = -1;
+        }
     }
 
-    
     public void process(
             EventHeader event) {
-        if (nevt == 0)
+        if (nevt == 0) {
             try {
 //open things up
                 fw = new FileWriter(outputTextName);
@@ -82,41 +101,49 @@
             } catch (IOException ex) {
                 Logger.getLogger(FastTrackAnalysisDriver.class.getName()).log(Level.SEVERE, null, ex);
             }
+        }
         //  Increment the event counter
         nevt++;
 
-      
 
-        Hep3Vector IP = new BasicHep3Vector(0., 0., 0.);
-        double bfield = event.getDetector().getFieldMap().getField(IP).z();
+
+        Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
+        double bfield = event.getDetector().getFieldMap().getField(IP).y();
 
         List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, "RawTrackerHitMaker_RawTrackerHits");
         List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
 
         List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
-        List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
+//        List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
 
-        int nAxialHitsTotal = axialhits.size();
+        int nAxialHitsTotal = stripHits.size();
         int nL1Hits = 0;
-        for (HelicalTrackHit hth : axialhits)
-            if (hth.Layer() == 1) nL1Hits++;
+        for (SiTrackerHitStrip1D str : stripHits) {
+            if (str.getRawHits().get(0).getLayerNumber()== 1) {
+                nL1Hits++;
+            }
+        }
 
 
         //  Create a relational table that maps TrackerHits to MCParticles
         RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
         List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
 
-        for (LCRelation relation : mcrelations)
-            if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+        for (LCRelation relation : mcrelations) {
+            if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
                 hittomc.add(relation.getFrom(), relation.getTo());
-
+            }
+        }
+/*
         RelationalTable hittomcAxial = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
 //        List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
         List<LCRelation> mcrelationsAxial = event.get(LCRelation.class, "AxialTrackMCRelations");
-        for (LCRelation relation : mcrelationsAxial)
-            if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+        for (LCRelation relation : mcrelationsAxial) {
+            if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
                 hittomcAxial.add(relation.getFrom(), relation.getTo());
-
+            }
+        }
+*/
         //  Instantiate the class that determines if a track is "findable"
         FindableTrack findable = new FindableTrack(event);
 
@@ -134,10 +161,13 @@
         List<List<SimTrackerHit>> simcols = event.get(SimTrackerHit.class);
 
         //  Loop over the SimTrackerHits and fill in the relational table
-        for (List<SimTrackerHit> simlist : simcols)
-            for (SimTrackerHit simhit : simlist)
-                if (simhit.getMCParticle() != null)
+        for (List<SimTrackerHit> simlist : simcols) {
+            for (SimTrackerHit simhit : simlist) {
+                if (simhit.getMCParticle() != null) {
                     mcHittomcP.add(simhit, simhit.getMCParticle());
+                }
+            }
+        }
 
         Map<Track, TrackAnalysis> tkanalMap = new HashMap<Track, TrackAnalysis>();
         Map<Track, BilliorTrack> btMap = new HashMap<Track, BilliorTrack>();
@@ -149,10 +179,12 @@
         int _nposRec = 0;
         //  Analyze the tracks in the event
         for (Track track : tracklist) {
-            if (track.getCharge() < 0)
+            if (track.getCharge() < 0) {
                 _neleRec++;
-            if (track.getCharge() > 0)
+            }
+            if (track.getCharge() > 0) {
                 _nposRec++;
+            }
             //extrapolate straight back...
             SeedTrack stEle = (SeedTrack) track;
             SeedCandidate seedEle = stEle.getSeedCandidate();
@@ -167,10 +199,11 @@
 
 
             //  Now analyze MC Particles on this track
-            MCParticle mcp = tkanal.getMCParticle();
-            if (mcp != null)
-                //  Create a map between the tracks found and the assigned MC particle
-                trktomc.add(track, tkanal.getMCParticle());
+            MCParticle mcp = tkanal.getMCParticleNew();
+            if (mcp != null) //  Create a map between the tracks found and the assigned MC particle
+            {
+                trktomc.add(track, tkanal.getMCParticleNew());
+            }
             BasicHep3Vector axial = new BasicHep3Vector();
             axial.setV(0, 1, 0);
             List<TrackerHit> hitsOnTrack = track.getTrackerHits();
@@ -192,20 +225,30 @@
                 double zerr = Math.sqrt(cross.getCorrectedCovMatrix().e(2, 2));
 
                 int htlayer = htc.Layer();
-                if (htlayer == 1) zlist[0] = z;
-                if (htlayer == 3) zlist[1] = z;
-                if (htlayer == 5) zlist[2] = z;
-                if (htlayer == 1) l1DeltaZ.put(track, z - zTr);
+                if (htlayer == 1) {
+                    zlist[0] = z;
+                }
+                if (htlayer == 3) {
+                    zlist[1] = z;
+                }
+                if (htlayer == 5) {
+                    zlist[2] = z;
+                }
+                if (htlayer == 1) {
+                    l1DeltaZ.put(track, z - zTr);
+                }
                 for (HelicalTrackStrip cl : clusterlist) {
                     int layer = cl.layer();
                     HelicalTrackStrip nearest = getNearestHit(cl, toththits);
                     if (layer == 1) {
                         Double l1Dist = getNearestDistance(cl, toththits);
-                        if (l1Dist != null)
+                        if (l1Dist != null) {
                             l1Isolation.put(track, l1Dist);
+                        }
                     }
-                    if (nearest != null)
+                    if (nearest != null) {
                         nearestHitToTrack.add(track, nearest);
+                    }
                 }
             }
             double layerDist = 10;//cm
@@ -229,9 +272,10 @@
             TrackAnalysis tkanal = tkanalMap.get(track);
             BilliorTrack bt = btMap.get(track);
             //  Calculate purity and make appropriate plots
-            MCParticle mcp = tkanal.getMCParticle();
-            if (mcp == null)
+            MCParticle mcp = tkanal.getMCParticleNew();
+            if (mcp == null) {
                 continue;
+            }
             if (mcp.getParents().size() == 1 && mcp.getParents().get(0).getPDGID() == 622) {
                 int nhits = tkanal.getNHitsNew();
                 double px = track.getPX();
@@ -260,11 +304,11 @@
         for (Track track1 : tracklist) {
             Track ele = null;
             Track pos = null;
-            int ch1 = track1.getCharge()*flipSign;
+            int ch1 = track1.getCharge() * flipSign;
             int index = tracklist.indexOf(track1);
             List<Track> subtracklist = tracklist.subList(index, tracklist.size());
             for (Track track2 : subtracklist) {
-                int ch2 = track2.getCharge()*flipSign;
+                int ch2 = track2.getCharge() * flipSign;
                 if (track1 != track2 && ch1 == -ch2) {
                     ele = track1;
                     pos = track2;
@@ -410,10 +454,12 @@
 
                     double l1minE = -99;
                     double l1minP = -99;
-                    if (l1Isolation.get(ele) != null)
+                    if (l1Isolation.get(ele) != null) {
                         l1minE = l1Isolation.get(ele);
-                    if (l1Isolation.get(pos) != null)
+                    }
+                    if (l1Isolation.get(pos) != null) {
                         l1minP = l1Isolation.get(pos);
+                    }
 
                     TrackAnalysis tkanalEle = tkanalMap.get(ele);
                     TrackAnalysis tkanalPos = tkanalMap.get(pos);
@@ -445,8 +491,12 @@
 
                     int eleFromAp = 0;
                     int posFromAp = 0;
-                    if (eleMC != null && ele == eleID) eleFromAp = 1;
-                    if (posMC != null && pos == posID) posFromAp = 1;
+                    if (eleMC != null && ele == eleID) {
+                        eleFromAp = 1;
+                    }
+                    if (posMC != null && pos == posID) {
+                        posFromAp = 1;
+                    }
                     MCParticle mcEle = tkanalEle.getMCParticle();
                     MCParticle mcPos = tkanalPos.getMCParticle();
                     double[] pmcEle = {-99, -99, -99};
@@ -472,18 +522,26 @@
                     //print out the vertex;
                     double l1dzE = -99;
                     double l1dzP = -99;
-                    if (l1DeltaZ.get(ele) != null)
+                    if (l1DeltaZ.get(ele) != null) {
                         l1dzE = l1DeltaZ.get(ele);
-                    if (l1DeltaZ.get(pos) != null)
+                    }
+                    if (l1DeltaZ.get(pos) != null) {
                         l1dzP = l1DeltaZ.get(pos);
+                    }
                     double l123KinkE = -99;
                     double l123KinkP = -99;
-                    if (l123KinkAngle.get(ele) != null)
+                    if (l123KinkAngle.get(ele) != null) {
                         l123KinkE = l123KinkAngle.get(ele);
-                    if (l123KinkAngle.get(pos) != null)
+                    }
+                    if (l123KinkAngle.get(pos) != null) {
                         l123KinkP = l123KinkAngle.get(pos);
-                    if (!tkanalEle.hasLayerOne()) nElectron = -nElectron;
-                    if (!tkanalPos.hasLayerOne()) nPositron = -nPositron;
+                    }
+                    if (!tkanalEle.hasLayerOne()) {
+                        nElectron = -nElectron;
+                    }
+                    if (!tkanalPos.hasLayerOne()) {
+                        nPositron = -nPositron;
+                    }
                     pw.format("%d %5.5f %5.5f %5.5f ", nevt, pxE, pyE, pzE);
                     pw.format("%5.5f  %5.5f %5.5f %5.5f  %5.5f ", d0E, z0E, slopeE, phi0E, RE);
                     pw.format("%5.5f  %5.5f %5.5f %5.5f ", chisqE, l1minE, l1dzE, l123KinkE);
@@ -561,46 +619,67 @@
         MCParticle mcApr = null;
         for (MCParticle mcp : mclist) {
 //            System.out.print("PDG ID = "+mcp.getPDGID());
-            if (mcp.getPDGID() == 622)
+            if (mcp.getPDGID() == 622) {
                 mcApr = mcp;
+            }
             if (mcp.getParents().size() == 1 && mcp.getParents().get(0).getPDGID() == 622) {
-                if (mcp.getPDGID() == -11) mcPos = mcp;
-                if (mcp.getPDGID() == 11) mcEle = mcp;
+                if (mcp.getPDGID() == -11) {
+                    mcPos = mcp;
+                }
+                if (mcp.getPDGID() == 11) {
+                    mcEle = mcp;
+                }
             }
         }
 //should probably check that each MC particles are found, but they should be...
-        double pxMCE=-99;
-        double pyMCE=-99;
-        double pzMCE=-99;
-        double pxMCP=-99;
-        double pyMCP=-99;
-        double pzMCP=-99;
+        double pxMCE = -99;
+        double pyMCE = -99;
+        double pzMCE = -99;
+        double pxMCP = -99;
+        double pyMCP = -99;
+        double pzMCP = -99;
+        double pxMCA = -99;
+        double pyMCA = -99;
+        double pzMCA = -99;
+        double mMCA = -99;
+        Hep3Vector decayMCA=new BasicHep3Vector();
         int findableE = 0;
         int findableP = 0;
-        int foundE=0;
-        int foundP=0;
+        int foundE = 0;
+        int foundP = 0;
+       if(mcEle!=null){
         pxMCE = mcEle.getPX();
-         pyMCE = mcEle.getPY();
-         pzMCE = mcEle.getPZ();
-        if (findable.InnerTrackerIsFindable(mcEle, _minLayers))
+        pyMCE = mcEle.getPY();
+        pzMCE = mcEle.getPZ();
+        if (findable.InnerTrackerIsFindable(mcEle, _minLayers)) {
             findableE = 1;
+        }
         Set<Track> trklistE = trktomc.allTo(mcEle);
         foundE = trklistE.size();//can be greater than 1 if more than 1 track shares hits
-
+       }else{
+           System.out.println("!!!!!   WHAT, no mcEle      !!!!!!!!!");
+       }
+       if(mcPos!=null){
         pxMCP = mcPos.getPX();
         pyMCP = mcPos.getPY();
         pzMCP = mcPos.getPZ();
-        if (findable.InnerTrackerIsFindable(mcPos, _minLayers))
+        if (findable.InnerTrackerIsFindable(mcPos, _minLayers)) {
             findableP = 1;
+        }
         Set<Track> trklistP = trktomc.allTo(mcPos);
         foundP = trklistP.size();//can be greater than 1 if more than 1 track shares hits
-
-        double pxMCA = mcApr.getPX();
-        double pyMCA = mcApr.getPY();
-        double pzMCA = mcApr.getPZ();
-        double mMCA = mcApr.getMass();
-        Hep3Vector decayMCA = mcApr.getEndPoint();
-
+       }else{
+            System.out.println("!!!!!   WHAT, no mcPos      !!!!!!!!!");
+       }
+       if(mcApr!=null){
+         pxMCA = mcApr.getPX();
+         pyMCA = mcApr.getPY();
+         pzMCA = mcApr.getPZ();
+        mMCA = mcApr.getMass();
+         decayMCA = mcApr.getEndPoint();
+       }else{
+           System.out.println("!!!!!   WHAT, no mvApr      !!!!!!!!!");
+       }
         pw.format(
                 "%d %5.5f %5.5f %5.5f %d %d ", -666, pxMCE, pyMCE, pzMCE, findableE, foundE);
         pw.format(
@@ -616,7 +695,6 @@
 
     }
 
-    
     public void endOfData() {
 
         pw.close();
@@ -652,7 +730,9 @@
         double dzdx1 = slt1.dzdx();
         double s1sq = 1 + 1 / (dydx1 * dydx1) + (dzdx1 * dzdx1) / (dydx1 * dydx1);
         double truep1y = Math.sqrt(p1mag2 / s1sq);
-        if (dydx1 < 0) truep1y = -truep1y;
+        if (dydx1 < 0) {
+            truep1y = -truep1y;
+        }
         double truep1x = truep1y / dydx1;
         double truep1z = dzdx1 * truep1x;
 
@@ -666,7 +746,9 @@
         double dzdx2 = slt2.dzdx();
         double s2sq = 1 + 1 / (dydx2 * dydx2) + (dzdx2 * dzdx2) / (dydx2 * dydx2);
         double truep2y = Math.sqrt(p2mag2 / s2sq);
-        if (dydx2 < 0) truep2y = -truep2y;
+        if (dydx2 < 0) {
+            truep2y = -truep2y;
+        }
         double truep2x = truep2y / dydx2;
         double truep2z = dzdx2 * truep2x;
 
@@ -701,7 +783,6 @@
         return doca;
     }
 
-
 //find the XOCA to the beamline extrpolating linearly from the reference point
     private double findXoca(double y, double z, double px, double py, double pz) {
         double xoca = 0;
@@ -754,7 +835,9 @@
         double dzdx1 = slt1.dzdx();
         double s1sq = 1 + 1 / (dydx1 * dydx1) + (dzdx1 * dzdx1) / (dydx1 * dydx1);
         truep[1] = Math.sqrt(p1mag2 / s1sq);
-        if (dydx1 < 0) truep[1] = -truep[1];
+        if (dydx1 < 0) {
+            truep[1] = -truep[1];
+        }
         truep[0] = truep[1] / dydx1;
         truep[2] = dzdx1 * truep[0];
         return new BasicHep3Vector(truep[0], truep[1], truep[2]);
@@ -779,7 +862,7 @@
         double mindist = 99999999;
         for (HelicalTrackHit hth : toththits) {
             HelicalTrackCross cross = (HelicalTrackCross) hth;
-            for (HelicalTrackStrip str : cross.getStrips())
+            for (HelicalTrackStrip str : cross.getStrips()) {
                 if (str != cl) {
                     Hep3Vector strorigin = str.origin();
                     Hep3Vector stru = str.u();
@@ -794,6 +877,7 @@
                     }
 
                 }
+            }
         }
 
         return nearest;
@@ -810,7 +894,7 @@
         Double mindist = 99999999.0;
         for (HelicalTrackHit hth : toththits) {
             HelicalTrackCross cross = (HelicalTrackCross) hth;
-            for (HelicalTrackStrip str : cross.getStrips())
+            for (HelicalTrackStrip str : cross.getStrips()) {
                 if (str != cl) {
                     Hep3Vector strorigin = str.origin();
                     Hep3Vector stru = str.u();
@@ -820,35 +904,40 @@
                     int strlayer = str.layer();
                     if (layer == strlayer && VecOp.sub(clvec, strvec).magnitude() < Math.abs(mindist)) {
                         mindist = VecOp.sub(clvec, strvec).magnitude();
-                        if (Math.abs(clvec.z()) > Math.abs(strvec.z()))
+                        if (Math.abs(clvec.z()) > Math.abs(strvec.z())) {
                             mindist = -mindist;
+                        }
                         nearest =
                                 str;
                     }
 
                 }
+            }
         }
         return mindist;
     }
 
     private Integer encodeBadHitList(List<Integer> badHits) {
         Integer badHitsEncoded = 0;
-        for (Integer layer : badHits)
+        for (Integer layer : badHits) {
             badHitsEncoded += (int) Math.pow(2, layer - 1);
+        }
         return badHitsEncoded;
     }
 
     private Integer encodeSharedHitList(List<Integer> sharedHits) {
         Integer sharedHitsEncoded = 0;
-        for (Integer layer : sharedHits)
+        for (Integer layer : sharedHits) {
             sharedHitsEncoded += (int) Math.pow(2, layer - 1);
+        }
         return sharedHitsEncoded;
     }
 
     private double getErr(double errSquared) {
-        if (errSquared > 0) return Math.sqrt(errSquared);
-        else return -99;
+        if (errSquared > 0) {
+            return Math.sqrt(errSquared);
+        } else {
+            return -99;
+        }
     }
 }
-
-

hps-java/src/main/java/org/lcsim/hps/users/mgraham
DetailedAnalysisDriver.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- DetailedAnalysisDriver.java	1 Aug 2011 18:25:20 -0000	1.2
+++ DetailedAnalysisDriver.java	29 Aug 2012 15:40:46 -0000	1.3
@@ -54,7 +54,7 @@
 
 /**
  *
- * @author partridge
+ * @author mgraham
  */
 public class DetailedAnalysisDriver extends Driver {
 
@@ -99,7 +99,7 @@
     public String outputPlots = "myplots.aida";
     Map<String, IProfile1D> clsizeMap = new HashMap<String, IProfile1D>();
     String[] detNames = {"Tracker"};
-    Integer[] nlayers = {8};
+    Integer[] nlayers = {10};
     int trk_count = 0;
     int nevt = 0;
     int _nmcTrk = 0;
@@ -141,14 +141,14 @@
         ctheffElectrons = hf.createProfile1D("Electrons Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
         d0effElectrons = hf.createProfile1D("Electrons Efficiency vs d0", "", 20, -1., 1.);
         z0effElectrons = hf.createProfile1D("Electrons Efficiency vs z0", "", 20, -1., 1.);
-
+/*
         peffAxial = hf.createProfile1D("Axial Efficiency vs p", "", 20, 0., beamP);
         thetaeffAxial = hf.createProfile1D("Axial Efficiency vs theta", "", 20, 80, 100);
         phieffAxial = hf.createProfile1D("Axial Efficiency vs phi", "", 25, -0.25, 0.25);
         ctheffAxial = hf.createProfile1D("Axial Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
         d0effAxial = hf.createProfile1D("Axial Efficiency vs d0", "", 20, -1., 1.);
         z0effAxial = hf.createProfile1D("Axial Efficiency vs z0", "", 20, -1., 1.);
-
+*/
         cthfake = hf.createProfile1D("Fake rate vs  cos(theta)", "", 25, -0.25, 0.25);
         phifake = hf.createProfile1D("Fake rate vs phi", "", 25, -0.25, 0.25);
         pfake = hf.createProfile1D("Fake rate vs p", "", 20, 0, 6);
@@ -202,8 +202,8 @@
         String debugDir = "debugPlots/";
         String occDir = "occupancyPlots/";
         //  Get the magnetic field
-        Hep3Vector IP = new BasicHep3Vector(0., 0., 0.);
-        double bfield = event.getDetector().getFieldMap().getField(IP).z();
+       Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
+        double bfield = event.getDetector().getFieldMap().getField(IP).y();
 
         List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, "RawTrackerHitMaker_RawTrackerHits");
         List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
@@ -279,14 +279,15 @@
 
 
 //        List<HelicalTrackHit> hthits = event.get(HelicalTrackHit.class, "MatchedHTHits");
-        List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
-        List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
+//        List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
+        List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "RotatedHelicalTrackHits");
+   //    List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
 
-        int nAxialHitsTotal = axialhits.size();
+//        int nAxialHitsTotal = axialhits.size();
         int nL1Hits = 0;
-        for (HelicalTrackHit hth : axialhits)
-            if (hth.Layer() == 1)
-                nL1Hits++;
+//        for (HelicalTrackHit hth : axialhits)
+//            if (hth.Layer() == 1)
+//                nL1Hits++;
 
         Map<String, Integer> occupancyMap = new HashMap<String, Integer>();
         for (RawTrackerHit rh : rawHits) {
@@ -332,28 +333,29 @@
 
         //  Create a relational table that maps TrackerHits to MCParticles
         RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
-        List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+       // List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+        List<LCRelation> mcrelations = event.get(LCRelation.class, "RotatedMCRelations");
 
         for (LCRelation relation : mcrelations)
             if (relation != null && relation.getFrom() != null && relation.getTo() != null)
                 hittomc.add(relation.getFrom(), relation.getTo());
 
-        RelationalTable hittomcAxial = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
-//        List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
-        List<LCRelation> mcrelationsAxial = event.get(LCRelation.class, "AxialTrackMCRelations");
-        for (LCRelation relation : mcrelationsAxial)
-            if (relation != null && relation.getFrom() != null && relation.getTo() != null)
-                hittomcAxial.add(relation.getFrom(), relation.getTo());
+//        RelationalTable hittomcAxial = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+////        List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+//        List<LCRelation> mcrelationsAxial = event.get(LCRelation.class, "AxialTrackMCRelations");
+//        for (LCRelation relation : mcrelationsAxial)
+//            if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+//                hittomcAxial.add(relation.getFrom(), relation.getTo());
 
         //  Instantiate the class that determines if a track is "findable"
         FindableTrack findable = new FindableTrack(event);
 
         //  Create a map between tracks and the associated MCParticle
         List<Track> tracklist = event.get(Track.class, "MatchedTracks");
-        List<Track> axialtracklist = event.get(Track.class, "AxialTracks");
+     //   List<Track> axialtracklist = event.get(Track.class, "AxialTracks");
 
         aida.cloud1D("Matched Tracks per Event").fill(tracklist.size());
-        aida.cloud1D("Axial Tracks per Event").fill(axialtracklist.size());
+ //       aida.cloud1D("Axial Tracks per Event").fill(axialtracklist.size());
         aida.cloud1D("HelicalTrackHits per Event").fill(toththits.size());
         RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
         RelationalTable trktomcAxial = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
@@ -382,12 +384,12 @@
         Map<Track, Double> l1Isolation = new HashMap<Track, Double>();
         Map<Track, Double> l1DeltaZ = new HashMap<Track, Double>();
         Map<Track, BilliorTrack> btMap = new HashMap<Track, BilliorTrack>();
-        String atrackdir = "TrackInfoAxial/";
-        //  Analyze the tracks in the event
-        for (Track atrack : axialtracklist) {
-            double apx = atrack.getPX();
-            aida.cloud1D(atrackdir + "pX").fill(apx);
-        }
+//        String atrackdir = "TrackInfoAxial/";
+//        //  Analyze the tracks in the event
+//        for (Track atrack : axialtracklist) {
+//            double apx = atrack.getPX();
+//            aida.cloud1D(atrackdir + "pX").fill(apx);
+//        }
         String trackdir = "TrackInfo/";
         //  Analyze the tracks in the event
         for (Track track : tracklist) {
@@ -442,9 +444,9 @@
             int nhits = tkanal.getNHitsNew();
             double purity = tkanal.getPurityNew();
             int nbad = tkanal.getNBadHitsNew();
-            int nbadAxial = tkanal.getNBadAxialHits();
+          //  int nbadAxial = tkanal.getNBadAxialHits();
             int nbadZ = tkanal.getNBadZHits();
-            int nAxial = tkanal.getNAxialHits();
+           // int nAxial = tkanal.getNAxialHits();
             int nZ = tkanal.getNZHits();
             List<Integer> badLayers = tkanal.getBadHitList();
             Integer badLayerEle = encodeBadHitList(badLayers);
@@ -454,9 +456,9 @@
             }
             aida.cloud1D(trackdir + "Mis-matched hits for all tracks").fill(nbad);
             aida.cloud1D(trackdir + "purityNew for all tracks").fill(purity);
-            aida.cloud1D(trackdir + "Bad Axial hits for all tracks").fill(nbadAxial);
+           // aida.cloud1D(trackdir + "Bad Axial hits for all tracks").fill(nbadAxial);
             aida.cloud1D(trackdir + "Bad Z hits for all tracks").fill(nbadZ);
-            aida.cloud1D(trackdir + "Number of Axial hits for all tracks").fill(nAxial);
+          //  aida.cloud1D(trackdir + "Number of Axial hits for all tracks").fill(nAxial);
             aida.cloud1D(trackdir + "Number of Z hits for all tracks").fill(nZ);
 
             for (Integer bhit : badLayers)
@@ -499,11 +501,11 @@
 
 
             //  Now analyze MC Particles on this track
-            MCParticle mcp = tkanal.getMCParticle();
+            MCParticle mcp = tkanal.getMCParticleNew();
             if (mcp != null) {
 
                 //  Create a map between the tracks found and the assigned MC particle
-                trktomc.add(track, tkanal.getMCParticle());
+                trktomc.add(track, tkanal.getMCParticleNew());
 
                 //  Calculate the MC momentum and polar angle
                 Hep3Vector Pmc = mcp.getMomentum();
@@ -705,14 +707,14 @@
         }
 
 //analyze the axial only tracks
-        for (Track track : axialtracklist) {
-            TrackAnalysis tkanal = new TrackAnalysis(track, hittomcAxial);
-            MCParticle mcp = tkanal.getMCParticle();
-            if (mcp != null)
-                //  Create a map between the tracks found and the assigned MC particle
-                trktomcAxial.add(track, tkanal.getMCParticle());
-
-        }
+//        for (Track track : axialtracklist) {
+//            TrackAnalysis tkanal = new TrackAnalysis(track, hittomcAxial);
+//            MCParticle mcp = tkanal.getMCParticle();
+//            if (mcp != null)
+//                //  Create a map between the tracks found and the assigned MC particle
+//                trktomcAxial.add(track, tkanal.getMCParticle());
+//
+//        }
         for (HelicalTrackHit hit : toththits) {
 
             int nAssHits = hit.getRawHits().size();
@@ -803,7 +805,7 @@
 
             TrackAnalysis tkanal = tkanalMap.get(track);          
             //  Calculate purity and make appropriate plots
-            MCParticle mcp = tkanal.getMCParticle();
+            MCParticle mcp = tkanal.getMCParticleNew();
             if (mcp == null)
                 continue;
             if (mcp.getParents().size() == 1 && mcp.getParents().get(0).getPDGID() == 622) {
@@ -956,7 +958,8 @@
             //  Find the number of layers hit by this mc particle
 //            System.out.println("MC pt=" + pt);
             int nhits = findable.LayersHit(mcp);
-            boolean isFindable = findable.InnerTrackerIsFindable(mcp, nlayers[0] - 2);
+//            boolean isFindable = findable.InnerTrackerIsFindable(mcp, nlayers[0] - 2);
+               boolean isFindable = findable.InnerTrackerIsFindable(mcp, nlayers[0]);
             Set<SimTrackerHit> mchitlist = mcHittomcP.allTo(mcp);
 
             Set<HelicalTrackCross> hitlist = hittomc.allTo(mcp);
@@ -972,8 +975,8 @@
             int ntrk = trklist.size();
 
 
-            Set<Track> trklistAxial = trktomcAxial.allTo(mcp);
-            int ntrkAxial = trklistAxial.size();
+//            Set<Track> trklistAxial = trktomcAxial.allTo(mcp);
+//            int ntrkAxial = trklistAxial.size();
 
             if (mcp.getPDGID() == 622) {
                 boolean bothreco = true;
@@ -984,7 +987,8 @@
                     if (trktomc.allTo(d).isEmpty())
                         bothreco = false;
 //                    if (findable.LayersHit(d) != nlayers[0])
-                    if (!findable.InnerTrackerIsFindable(d, nlayers[0] - 2))
+//                    if (!findable.InnerTrackerIsFindable(d, nlayers[0] - 2))
+                     if (!findable.InnerTrackerIsFindable(d, nlayers[0]))
                         bothfindable = false;
                 }
                 double vtxWgt = 0;
@@ -1030,15 +1034,15 @@
                 //              if (wgt == 0)
 //                    System.out.println("Missed a findable track!");
 
-                double wgtAxial = 0.;
-                if (ntrkAxial > 0)
-                    wgtAxial = 1.;
-                peffAxial.fill(p, wgtAxial);
-                phieffAxial.fill(phi, wgtAxial);
-                thetaeffAxial.fill(theta, wgtAxial);
-                ctheffAxial.fill(cth, wgtAxial);
-                d0effAxial.fill(d0, wgtAxial);
-                z0effAxial.fill(z0, wgtAxial);
+//                double wgtAxial = 0.;
+//                if (ntrkAxial > 0)
+//                    wgtAxial = 1.;
+//                peffAxial.fill(p, wgtAxial);
+//                phieffAxial.fill(phi, wgtAxial);
+//                thetaeffAxial.fill(theta, wgtAxial);
+//                ctheffAxial.fill(cth, wgtAxial);
+//                d0effAxial.fill(d0, wgtAxial);
+//                z0effAxial.fill(z0, wgtAxial);
 
 
             }
@@ -1062,15 +1066,15 @@
 //                        System.out.println("Missed a findable ELECTRON!!!!!");
 
 
-                    double wgtAxial = 0.;
-                    if (ntrkAxial > 0)
-                        wgtAxial = 1.;
-                    peffAxial.fill(p, wgtAxial);
-                    phieffAxial.fill(phi, wgtAxial);
-                    thetaeffAxial.fill(theta, wgtAxial);
-                    ctheffAxial.fill(cth, wgtAxial);
-                    d0effAxial.fill(d0, wgtAxial);
-                    z0effAxial.fill(z0, wgtAxial);
+//                    double wgtAxial = 0.;
+//                    if (ntrkAxial > 0)
+//                        wgtAxial = 1.;
+//                    peffAxial.fill(p, wgtAxial);
+//                    phieffAxial.fill(phi, wgtAxial);
+//                    thetaeffAxial.fill(theta, wgtAxial);
+//                    ctheffAxial.fill(cth, wgtAxial);
+//                    d0effAxial.fill(d0, wgtAxial);
+//                    z0effAxial.fill(z0, wgtAxial);
 
                 }
 

hps-java/src/main/java/org/lcsim/hps/recon/tracking
FindableTrack.java 1.1 -> 1.2
diff -u -r1.1 -r1.2
--- FindableTrack.java	1 Aug 2011 18:25:21 -0000	1.1
+++ FindableTrack.java	29 Aug 2012 15:40:46 -0000	1.2
@@ -71,7 +71,7 @@
         _nlayersTot=ntot;
         //  Get the magnetic field
         Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
-        _bfield = event.getDetector().getFieldMap().getField(IP).z();
+         _bfield = event.getDetector().getFieldMap().getField(IP).y();
 
         //  Instantiate the hit identifier class
         _ID = new HitIdentifier();

hps-java/src/main/resources/org/lcsim/hps/steering
HPSTrackingBunchedMC.lcsim added at 1.1
diff -N HPSTrackingBunchedMC.lcsim
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ HPSTrackingBunchedMC.lcsim	29 Aug 2012 15:40:46 -0000	1.1
@@ -0,0 +1,58 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance" 
+       xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+       
+    <inputFiles>
+        <file>${inputFile}</file>
+    </inputFiles>
+    
+    <control>
+        <numberOfEvents>-1</numberOfEvents>
+        <printInputFiles>true</printInputFiles>
+        <printDriversDetailed>true</printDriversDetailed>
+    </control>
+    
+    <execute>
+        <driver name="EventMarkerDriver"/>
+        <driver name="SVTSetupDriver"/>
+        <driver name="TrackerDigiDriver"/>
+        <driver name="TrackerReconDriver"/>
+        <driver name="AnalysisDriver"/>
+<!--        <driver name="LCIOWriter"/>  -->
+    </execute>
+    
+    <drivers>
+
+        <driver name="SVTSetupDriver"
+                type="org.lcsim.hps.recon.tracking.HPSSVTSensorSetup">   
+        </driver>
+
+        <driver name="EventMarkerDriver"
+                type="org.lcsim.job.EventMarkerDriver">
+            <eventInterval>1</eventInterval>
+        </driver>
+    
+        <driver name="TrackerDigiDriver"
+                type="org.lcsim.hps.recon.tracking.TrackerDigiDriver">
+            <debug>false</debug>
+        </driver>
+                
+        <driver name="TrackerReconDriver"
+                type="org.lcsim.hps.recon.tracking.TrackerReconDriver">
+            <debug>false</debug>
+            <strategyResource>/org/lcsim/hps/recon/tracking/strategies/HPS-Test-4pt0.xml</strategyResource>
+            <stripMaxSeparation>20.0</stripMaxSeparation>
+            <stripTolerance>1.0</stripTolerance>
+        </driver>
+
+        <driver name="AnalysisDriver" type="org.lcsim.hps.users.mgraham.FastTrackAnalysisDriver">
+            <outputText>${output}</outputText>
+        </driver>
+        <!-- 
+        <driver name="LCIOWriter"
+                type="org.lcsim.util.loop.LCIODriver">
+            <outputFilePath>${inputFile}_trackRecon</outputFilePath>
+        </driver>
+                          -->                                          
+    </drivers>
+</lcsim>

hps-java/src/main/java/org/lcsim/hps/recon/vertexing
BilliorVertex.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- BilliorVertex.java	7 Jul 2011 20:57:39 -0000	1.2
+++ BilliorVertex.java	29 Aug 2012 15:40:46 -0000	1.3
@@ -11,7 +11,7 @@
 import static java.lang.Math.tan;
 
 /**
- * @version $Id: BilliorVertex.java,v 1.2 2011/07/07 20:57:39 mgraham Exp $
+ * @version $Id: BilliorVertex.java,v 1.3 2012/08/29 15:40:46 mgraham Exp $
  * @version Vertex tracks using least-squares method laid out by billior etal used in the HPS Java package.
  */
 public class BilliorVertex {
@@ -260,7 +260,7 @@
         }
         _covVtx = MatrixOp.inverse(sumwi);
         if (_debug)
-            System.out.println("_covVtx matrix " + _covVtx.toString());
+            System.out.println("fastVertex::_covVtx matrix " + _covVtx.toString());
         _vertexPosition = (BasicMatrix) MatrixOp.mult(_covVtx, sumwiXi);
         _chiSq = 0;
         //get the chisq
@@ -359,7 +359,7 @@
             WList.add(tmpW);
 
             if (_debug)
-                System.out.println("Params = \n" + tmpPar);
+                System.out.println("makeDerivativeMatrices::Params = \n" + tmpPar);
             if (_debug)
                 System.out.println("D = \n" + tmpD);
             if (_debug)
@@ -388,7 +388,7 @@
         int n = 1;
         for (Matrix covVtxMom : covVtxMomList) {
             if (_debug)
-                System.out.println("Track " + n + "  covVtxMom : " + covVtxMom.toString());
+                System.out.println("addV0fromBSConstraint::Track " + n + "  covVtxMom : " + covVtxMom.toString());
             MatrixOp.setSubMatrix(Ckm1, covVtxMom, 0, 3 * n);
             MatrixOp.setSubMatrix(Ckm1, MatrixOp.transposed(covVtxMom), 3 * n, 0);
             n++;
@@ -397,7 +397,7 @@
             BasicMatrix pi = (BasicMatrix) _pFit.get(i);
             MatrixOp.setSubMatrix(Xkm1, pi, 3 * (i + 1), 0);
             if (_debug)
-                System.out.println("Track " + i + "  p : " + pi.toString());
+                System.out.println("addV0fromBSConstraint::Track " + i + "  p : " + pi.toString());
             for (int j = 0; j < _ntracks; j++)
                 MatrixOp.setSubMatrix(Ckm1, covMomList[i][j], 3 * (i + 1), 3 * (j + 1));
         }
@@ -430,7 +430,7 @@
         //calculate the position of the A' at X=0
         BasicMatrix rk = new BasicMatrix(3, 1);
         if (_debug)
-            System.out.println("Vx = " + Vx + "; Vy = " + Vy + "; Vz = " + Vz + "; pxtot = " + pxtot + "; pytot = " + pytot + "; pztot = " + pztot);
+            System.out.println("addV0fromBSConstraint::Vx = " + Vx + "; Vy = " + Vy + "; Vz = " + Vz + "; pxtot = " + pxtot + "; pytot = " + pytot + "; pztot = " + pztot);
         rk.setElement(0, 0, 0);
         rk.setElement(1, 0, 0 - (Vy - pytot / pxtot * Vx));
         rk.setElement(2, 0, 0 - (Vz - pztot / pxtot * Vx));
@@ -487,7 +487,7 @@
         //now do the matrix operations to get the constrained parameters
         BasicMatrix Hkt = (BasicMatrix) MatrixOp.transposed(Hk);
         if (_debug)
-            System.out.println("Ckm1Hk = " + MatrixOp.mult(Ckm1, Hk));
+            System.out.println("addV0fromBSConstraint::Ckm1Hk = " + MatrixOp.mult(Ckm1, Hk));
 
         BasicMatrix Rk = (BasicMatrix) MatrixOp.mult(Hkt, MatrixOp.mult(Ckm1, Hk));
         if (_debug)
@@ -548,7 +548,7 @@
         int n = 1;
         for (Matrix covVtxMom : covVtxMomList) {
             if (_debug)
-                System.out.println("Track " + n + "  covVtxMom : " + covVtxMom.toString());
+                System.out.println("constrainV0toBS::Track " + n + "  covVtxMom : " + covVtxMom.toString());
             MatrixOp.setSubMatrix(Ckm1, covVtxMom, 0, 3 * n);
             MatrixOp.setSubMatrix(Ckm1, MatrixOp.transposed(covVtxMom), 3 * n, 0);
             n++;
@@ -732,9 +732,12 @@
         mom[0] = Pt * Math.cos(phiv);
         mom[1] = Pt * Math.sin(phiv);
         mom[2] = Pt * 1 / Math.tan(theta);
-        if (_debug)
-            System.out.println(mom[0] + "; " + mom[1] + "; " + mom[2]);
-        return mom;
+        if (_debug){
+            System.out.println("getFittedMomentum::  "+mom[0] + "; " + mom[1] + "; " + mom[2]);
+
+            System.out.println("pT= "+Pt+"; phi = "+phiv+"; B = "+ _bField);
+        }
+      return mom;
     }
 
     public double getInvMass() {
@@ -899,7 +902,7 @@
         BasicMatrix covVtx = (BasicMatrix) MatrixOp.inverse(tmpCovVtx);
         BasicMatrix xtilde = (BasicMatrix) MatrixOp.mult(covVtx, bigsum);
         if (_debug)
-            System.out.println("Vertex at : \nx= " + xtilde.e(0, 0) + " +/- " + Math.sqrt(covVtx.e(0, 0)) + "\ny= " + xtilde.e(1, 0) + " +/- " + Math.sqrt(covVtx.e(1, 1)) + "\nz= " + xtilde.e(2, 0) + " +/- " + Math.sqrt(covVtx.e(2, 2)));
+            System.out.println("follow1985Paper::Vertex at : \nx= " + xtilde.e(0, 0) + " +/- " + Math.sqrt(covVtx.e(0, 0)) + "\ny= " + xtilde.e(1, 0) + " +/- " + Math.sqrt(covVtx.e(1, 1)) + "\nz= " + xtilde.e(2, 0) + " +/- " + Math.sqrt(covVtx.e(2, 2)));
 
         //ok, now the momentum
         List<Matrix> qtildes = new ArrayList<Matrix>();
@@ -926,7 +929,7 @@
             ptildes.add(ptilde);
             chisq += MatrixOp.mult(MatrixOp.transposed(MatrixOp.add(p, MatrixOp.mult(-1, ptilde))), MatrixOp.mult(g, MatrixOp.add(p, MatrixOp.mult(-1, ptilde)))).e(0, 0);
             if (_debug)
-                System.out.println("\n\nTrack #" + j);
+                System.out.println("\n\nfollow1985Paper::Track #" + j);
             if (_debug)
                 System.out.println("eps(meas)    = " + p.e(0, 0) + "        eps(fit)   =" + ptilde.e(0, 0));
             if (_debug)
@@ -958,7 +961,7 @@
         }
 
         if (_debug)
-            System.out.println("chi^2 = " + chisq);
+            System.out.println("follow1985Paper::chi^2 = " + chisq);
 
         _chiSq = chisq;
         _covVtx = covVtx;

hps-java/src/main/java/org/lcsim/hps/util
ConvertToStdhep.java 1.4 -> 1.5
diff -u -r1.4 -r1.5
--- ConvertToStdhep.java	6 Mar 2012 05:12:15 -0000	1.4
+++ ConvertToStdhep.java	29 Aug 2012 15:40:46 -0000	1.5
@@ -18,10 +18,21 @@
 import hep.io.stdhep.StdhepEndRun;
 import hep.io.stdhep.StdhepEvent;
 import hep.io.stdhep.StdhepWriter;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
 import java.io.File;
 import java.util.Random;
 import java.util.logging.Level;
 import java.util.logging.Logger;
+import org.apache.commons.cli.CommandLine;
+import org.apache.commons.cli.CommandLineParser;
+import org.apache.commons.cli.HelpFormatter;
+import org.apache.commons.cli.Option;
+import org.apache.commons.cli.Options;
+import org.apache.commons.cli.ParseException;
+import org.apache.commons.cli.PosixParser;
+import org.lcsim.detector.IRotation3D;
+import org.lcsim.detector.RotationGeant;
 
 public class ConvertToStdhep {
 
@@ -32,58 +43,133 @@
     static double _xoff = 0.0;   //set the x,y,z origin offsets...
     static double _yoff = 0.0;
     static double _zoff = 0.03;
-    static double sigx = 0.01;        // Luminous region size in x (mm)
-    static double sigy = 0.01;        // Luminous region size in y (mm)
+    static double aMass = 99;  //Aprime mass (MeV)
+    static double sigx = 0.02;        // Luminous region size in x (mm)
+    static double sigy = 0.02;        // Luminous region size in y (mm)
     static double sigz = 0.0;         // Luminous region size in z (mm)
-    static double aMass = 0.05;  //Aprime mass (GeV)
+//beam is positioned so that at first beam direction is in z, then rotated to correct orientation    
+    static double rotx = 0.00;        // Rotation of beam about X
+    static double roty = 0.00;        // Rotation of beam about Y
+    static double rotz = 0.00;         // Rotation of beam about Z
     static double maxLen = 200; // maximum decay length (mm)
-    static double _ecm = 5.5; //GeV
-    static int nInFiles = 100;  //number of input files
+    static double _ecm = 2.2; //GeV
+    static int nInFiles = 10000;  //number of input files
     static StdhepWriter sw;
     static int nmax = 500000;  //maximum number of events to write to 1 stdhep file (new one opens if n>nmax)
     static String fileType = "lhe";
+    static IRotation3D rot = new RotationGeant(rotx, roty, rotz);
 //        static String fileType="dat";
+    static int nread = -1;  // a running total of number of events read/written to stdhep files
+    static boolean _eventFilter=false;
+    static int _nEleRequired=2;
+    
+    private static Options createCommandLineOptions() {
+        Options options = new Options();
+
+        options.addOption(new Option("m", true, "A' Mass (MeV)"));
+        options.addOption(new Option("e", true, "Beam Energy (GeV)"));
+        options.addOption(new Option("n", true, "Number of files to run."));
+         options.addOption(new Option("s", false, "Filter Events"));
+
+        return options;
+    }
 
     /**
      * @param args the command line arguments
      * @throws IOException
      */
     public static void main(String[] args) throws IOException {
-       
+
+
+        // Set up command line parsing.
+        Options options = createCommandLineOptions();
+
+        CommandLineParser parser = new PosixParser();
+
+        // Parse command line arguments.
+        CommandLine cl = null;
+        try {
+            cl = parser.parse(options, args);
+            System.out.println("Trying parser");
+        } catch (ParseException e) {
+            throw new RuntimeException("Problem parsing command line options.", e);
+        }
+        String massString = String.valueOf(aMass);
+        String ninString = String.valueOf(nInFiles);
+        String ecmString = String.valueOf(_ecm);
+        String eptString = String.valueOf(_ecm);
+        eptString = convertDecimal(eptString);
+        // LCIO output file.
+        if (cl.hasOption("m")) {
+            massString = cl.getOptionValue("m");
+            aMass = Integer.valueOf(massString);
+            System.out.println(massString);
+        }
+
+        if (cl.hasOption("n")) {
+            ninString = cl.getOptionValue("n");
+            nInFiles = Integer.valueOf(ninString);
+            System.out.println(ninString);
+        }
+
+        if (cl.hasOption("e")) {
+            ecmString = cl.getOptionValue("e");
+            _ecm = Double.valueOf(ninString);
+            System.out.println(ecmString);
+            eptString = convertDecimal(ecmString);
+        }
         
-        String postfix = "_20ux200u_beamspot_gammactau_0cm.stdhep";
+        String filter="all";
+        if (cl.hasOption("s")) {
+           _eventFilter=true;
+           filter="selected";
+        }
 
+        String postfix = "_20ux200u_beamspot_gammactau_0cm.stdhep";
 //        String fDir="/nfs/slac/g/hps/mgraham/DarkPhoton/tvm/testrun/";
 //        String fileLabel = "ap2.2gev40mevsel";
 
-  //            String fDir="/nfs/slac/g/hps/mgraham/DarkPhoton/tvm/trident/full/";
-  //      String fileLabel = "full6.6gev";
+        //            String fDir="/nfs/slac/g/hps/mgraham/DarkPhoton/tvm/trident/full/";
+        //      String fileLabel = "full6.6gev";
+
+        //String fDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/MadGraph/aMassEvents2pt2Ap100MeV/";
+        //String fileLabel = "ap2.2gev100mevall";
+        //String inLabel = "W2pt2GeV_Ap100MeV_";
+
+        String fDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/MadGraph/Events" + eptString + "Ap" + massString + "MeV/";
+        String fileLabel = "ap" + ecmString + "gev" + massString + "mev"+filter;
+        String inLabel = "W" + eptString + "GeV_Ap" + massString + "MeV_";
 
-       String fDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/MadGraph/Events1pt1Ap20MeV/";
-        String fileLabel = "ap1.1gev20mevall";
-        String inLabel = "W1pt1GeV_Ap20MeV_";
         String inPost = "_unweighted_events.lhe";
 
 
         String outDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/SignalEvents/";
-        int nread = -1;
+
         int nOutFile = 0;
-        
+        System.out.println(fDir + fileLabel);
         for (int i = 0; i < nInFiles; i++) {
             int fnum = i + 1;
             String snum = "_" + fnum;
-             if (fnum < 10000) snum = "_0" + fnum;
-            if (fnum < 1000) snum = "_00" + fnum;
-            if (fnum < 100) snum = "_000" + fnum;
-            if (fnum < 10) snum = "_0000" + fnum;
+            if (fnum < 10000) {
+                snum = "_0" + fnum;
+            }
+            if (fnum < 1000) {
+                snum = "_00" + fnum;
+            }
+            if (fnum < 100) {
+                snum = "_000" + fnum;
+            }
+            if (fnum < 10) {
+                snum = "_0000" + fnum;
+            }
 
             String infile = "";
-            if (fileType.contains("dat")){
+            if (fileType.contains("dat")) {
                 infile = fDir + fileLabel + snum + ".dat";
 //                System.out.println(infile);
-            }
-            else if (fileType.contains("lhe")) {
+            } else if (fileType.contains("lhe")) {
                 infile = fDir + inLabel + i + inPost;
+                System.out.println("Unzipping " + infile);
                 String cmd = "gunzip " + infile + ".gz";
                 Process p = Runtime.getRuntime().exec(cmd);
                 try {
@@ -92,14 +178,16 @@
                     Logger.getLogger(ConvertToStdhep.class.getName()).log(Level.SEVERE, null, ex);
                 }
             }
-             File f = new File(infile);
-            if (nread == -1&&f.exists()) {
+            File f = new File(infile);
+            if (nread == -1 && f.exists()) {
                 nOutFile++;
                 String outfile = outDir + fileLabel + "_" + nOutFile + postfix;//replace .txt by .stdhep
                 System.out.println("==== processing " + infile + " into " + outfile + "  aP mass = " + aMass + "====");
                 openStdHepFile(outfile);
-            }           
-            if (f.exists()) nread += process(infile);
+            }
+            if (f.exists()) {
+                nread += process(infile);
+            }
             if (fileType.contains("lhe")) {
                 String cmd = "gzip " + infile;
                 Process p = Runtime.getRuntime().exec(cmd);
@@ -121,10 +209,12 @@
     private static int lineCounter(StreamTokenizer tok) throws IOException {
         int lines = 0;
         while (tok.nextToken() != tok.TT_EOF) {
-            if (tok.ttype == tok.TT_EOL)
+            if (tok.ttype == tok.TT_EOL) {
                 lines++;
-            if (tok.ttype == tok.TT_WORD && tok.sval.startsWith("nev"))
+            }
+            if (tok.ttype == tok.TT_WORD && tok.sval.startsWith("nev")) {
                 return lines;
+            }
         }
         //shouldn't get here...but maybe
         return lines;
@@ -140,12 +230,13 @@
     }
 
     private static int getnevts(StreamTokenizer lctok) throws IOException {
-    int nevts=-1;
-        if (fileType.contains("dat"))
+        int nevts = -1;
+        if (fileType.contains("dat")) {
             return lineCounter(lctok);
-        else if (fileType.contains("lhe")) {
-            while (nevts == -1)
+        } else if (fileType.contains("lhe")) {
+            while (nevts == -1) {
                 nevts = getNumberOfEvents(lctok);
+            }
             return nevts;
         }
         return nevts;
@@ -185,8 +276,8 @@
         float xsec = (float) 99999997952.;
         double rn1 = 12345321;
         double rn2 = 66666666;
-        StdhepBeginRun sb = new StdhepBeginRun(nreq, ngen, nwrit, ecm, xsec, rn1, rn2);
-        sw.writeRecord(sb);
+//        StdhepBeginRun sb = new StdhepBeginRun(nreq, ngen, nwrit, ecm, xsec, rn1, rn2);
+//        sw.writeRecord(sb);
 
 
         tok.resetSyntax();
@@ -200,48 +291,60 @@
         tok.whitespaceChars(0, ' ');
         tok.eolIsSignificant(true);
 
-      
+
 
         double[] beam = {0, 0, 0, 0};
         int nevhep = 0;
         for (int icross = 0; icross < nwrit; icross++) {
-            beam[0] = sigx * generator.nextGaussian() + _xoff;
-            beam[1] = sigy * generator.nextGaussian() + _yoff;
-            beam[2] = sigz * generator.nextGaussian() + _zoff;
+            Hep3Vector beamVecOrig =
+                    new BasicHep3Vector(sigx * generator.nextGaussian() + _xoff,
+                    sigy * generator.nextGaussian() + _yoff,
+                    sigz * generator.nextGaussian() + _zoff);
+            Hep3Vector beamVec = rot.rotated(beamVecOrig);
+            beam[0] = beamVec.x();
+            beam[1] = beamVec.y();
+            beam[2] = beamVec.z();
+
             double tmpDecLen = 0;
-           
-            if (fileType.contains("lhe"))
-                writeLHEEvent(tok, beam,  icross);
-            else if (fileType.contains("dat"))
+
+            if (fileType.contains("lhe")) {
+                writeLHEEvent(tok, beam, icross);
+            } else if (fileType.contains("dat")) {
                 writeDATEvent(tok, beam, tmpDecLen, icross);
+            }
 
         }
         fr.close();
-        int nreqe = (int) nevts;
-        int ngene = (int) nevts;
-        int nwrite = (int) nevts;
-        float ecme = (float) ecm;
-        float xsece = (float) 99999997952.;
-        double rn1e = 12345321;
-        double rn2e = 66666666;
-        StdhepEndRun se = new StdhepEndRun(nreqe, ngene, nwrite, ecme, xsece, rn1e, rn2e);
-        sw.writeRecord(se);
+
+
         return nwrit;
 
 
     }
 
     private static void closeStdHepFile() throws IOException {
-        if(sw!=null) sw.close();
-        System.out.println("Ok...done!");
+        if (sw != null) {
+            int nreqe = (int) nread;
+            int ngene = (int) nread;
+            int nwrite = (int) nread;
+            float ecme = (float) _ecm;
+            float xsece = (float) 99999997952.;
+            double rn1e = 12345321;
+            double rn2e = 66666666;
+            StdhepEndRun se = new StdhepEndRun(nreqe, ngene, nwrite, ecme, xsece, rn1e, rn2e);
+            sw.writeRecord(se);
+            sw.close();
+        }
+        System.out.println("Ok...done with closing!");
     }
 
     @SuppressWarnings("static-access")
     private static List<Double> getNumbersInLine(StreamTokenizer tok) throws IOException {
         List<Double> nums = new ArrayList<Double>();
         while (tok.nextToken() != tok.TT_EOF) {
-            if (tok.ttype == tok.TT_EOL)
+            if (tok.ttype == tok.TT_EOL) {
                 break;
+            }
             String tokVal = tok.sval;
 //            System.out.println(tokVal);
             nums.add(Double.valueOf(tokVal).doubleValue());
@@ -258,38 +361,43 @@
         int evts = -1;
         while (tok.nextToken() != tok.TT_EOF) {
 //            System.out.println(tok.toString());
-            if (tok.ttype == tok.TT_EOL)
+            if (tok.ttype == tok.TT_EOL) {
                 break;
+            }
             if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("Number")) {
-  //              System.out.println(tok.toString());
+                //              System.out.println(tok.toString());
                 fndNumber = true;
             }
-            if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("of"))
+            if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("of")) {
                 fndOf = true;
-            if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("Events"))
+            }
+            if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("Events")) {
                 fndEvents = true;
-            if (tok.ttype == tok.TT_NUMBER && fndEvents && fndOf && fndNumber)
+            }
+            if (tok.ttype == tok.TT_NUMBER && fndEvents && fndOf && fndNumber) {
                 evts = (int) tok.nval;
+            }
         }
         return evts;
     }
 
     @SuppressWarnings("static-access")
     private static void getToNextEvent(StreamTokenizer tok) throws IOException {
-        while (tok.nextToken() != tok.TT_EOF)
-            //            System.out.println(tok.toString());
+        while (tok.nextToken() != tok.TT_EOF) //            System.out.println(tok.toString());
+        {
             if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("<event>")) {
                 tok.nextToken();//get to the EOL
                 return;
             }
+        }
     }
 
     static private double expWeight(double x) {
         return Math.exp(-x / _declength);
     }
-    
+
     static private double expWeight(double x, double gamma) {
-        return Math.exp(-x / (gamma*_declength));
+        return Math.exp(-x / (gamma * _declength));
     }
 
     static private double findMaxWeight() {
@@ -300,7 +408,9 @@
         for (int i = 0; i < ntrials; i++) {
             double x = generator.nextDouble() * maxlength;
             double wght = expWeight(x);
-            if (wght > maxWeight) maxWeight = wght;
+            if (wght > maxWeight) {
+                maxWeight = wght;
+            }
         }
 
         return maxWeight;
@@ -314,11 +424,11 @@
         double tmpwght = 0;
         while (tmpwght < draw) {
             dl = generator.nextDouble() * maxlength;
-            tmpwght = expWeight(dl,gamma) / MaxWeight;
+            tmpwght = expWeight(dl, gamma) / MaxWeight;
         }
         return dl;
     }
-    
+
     static private double getDecayLength(double MaxWeight) {
         Random generator = new Random();
         double maxlength = maxLen;
@@ -336,8 +446,9 @@
         List<Double> vals = getNumbersInLine(tok);
 //            System.out.println(nums.toString());
 
-        if (vals.size() != 16)
+        if (vals.size() != 16) {
             throw new RuntimeException("Unexpected entry for event:  size = " + vals.size());
+        }
 
 //            int nhep = nums.get(0).intValue();
         int nhep = 4;  // two daughters of the A' and the recoil e-
@@ -358,13 +469,21 @@
         for (int npart = 0; npart < nhep; npart++) {
 
             isthep[i] = 1;
-            if (npart == 0) isthep[i] = 2;
+            if (npart == 0) {
+                isthep[i] = 2;
+            }
 //                 if (npart == 0) isthep[i] = 0;
 //                 if (npart == 0) isthep[i] = 1;
             idhep[i] = 622;
-            if (npart == 1) idhep[i] = -11;
-            if (npart == 2) idhep[i] = 11;
-            if (npart == 3) idhep[i] = 11;
+            if (npart == 1) {
+                idhep[i] = -11;
+            }
+            if (npart == 2) {
+                idhep[i] = 11;
+            }
+            if (npart == 3) {
+                idhep[i] = 11;
+            }
             jmohep[2 * i] = 0;
             jmohep[2 * i + 1] = 0;
             jdahep[2 * i] = 2;
@@ -381,49 +500,52 @@
                 jdahep[2 * i] = 0;
                 jdahep[2 * i + 1] = 0;
             }
-            for (int j = 0; j < 4; j++)
+            for (int j = 0; j < 4; j++) {
                 phep[5 * i + j] = vals.get(4 * i + j);
+            }
             phep[5 * i + 4] = 0.51109989000E-03;
             if (npart == 0 && !trident) {
-                phep[5 * i + 4] = aMass;
+                phep[5 * i + 4] = aMass / 1000.0;
                 /*
-                ApMom[0] = phep[5 * i + 2];
-                ApMom[1] = phep[5 * i + 1];
-                ApMom[2] = phep[5 * i + 0];
+                 ApMom[0] = phep[5 * i + 2];
+                 ApMom[1] = phep[5 * i + 1];
+                 ApMom[2] = phep[5 * i + 0];
                  * 
                  */
-                 ApMom[0] = phep[5 * i + 0];
+                ApMom[0] = phep[5 * i + 0];
                 ApMom[1] = phep[5 * i + 1];
                 ApMom[2] = phep[5 * i + 2];
             }
-            for (int j = 0; j < 4; j++)
+            for (int j = 0; j < 4; j++) {
                 vhep[4 * i + j] = beam[j];
+            }
             //decay the A' and daughters daughters at _declength
-            if (!trident && (npart == 1 || npart == 2 || npart == 0))
-                if (!expDecay && !flatDecay)
+            if (!trident && (npart == 1 || npart == 2 || npart == 0)) {
+                if (!expDecay && !flatDecay) {
                     vhep[4 * i + 0] = beam[0] + _declength;
-                else {
-                    double totApMom = Math.sqrt(ApMom[0] * ApMom[0] + ApMom[1] * ApMom[1] + ApMom[2] * ApMom[2]);                  
+                } else {
+                    double totApMom = Math.sqrt(ApMom[0] * ApMom[0] + ApMom[1] * ApMom[1] + ApMom[2] * ApMom[2]);
 //                        System.out.println("Decay at : " + tmpDecLen);
                     vhep[4 * i + 0] = beam[0] + decLen * ApMom[0] / totApMom;
                     vhep[4 * i + 1] = beam[1] + decLen * ApMom[1] / totApMom;
                     vhep[4 * i + 2] = beam[2] + decLen * ApMom[2] / totApMom;
                 }
+            }
             double px = phep[5 * i + 0];
             double pz = phep[5 * i + 2];
 //            phep[5 * i + 0] = pz;// these two are from back in the lcsim/jlab coordinate system days
 //            phep[5 * i + 2] = px;//
-             phep[5 * i + 0] = px;
+            phep[5 * i + 0] = px;
             phep[5 * i + 2] = pz;
             //                                if (i == 0 || i == nhep - 1) {
 /*
-            System.out.println(i + " st: " + isthep[i] + " id: " + idhep[i] +
-            " jmo: " + jmohep[2 * i] + " " + jmohep[2 * i + 1] +
-            " jda: " + jdahep[2 * i] + " " + jdahep[2 * i + 1]);
-            System.out.println("p: " + phep[5 * i] + " " + phep[5 * i + 1] + " " +
-            phep[5 * i + 2] + " " + phep[5 * i + 3] + " " + phep[5 * i + 4]);
-            System.out.println("v: " + vhep[4 * i] + " " + vhep[4 * i + 1] + " " +
-            vhep[4 * i + 2] + " " + vhep[4 * i + 3]);
+             System.out.println(i + " st: " + isthep[i] + " id: " + idhep[i] +
+             " jmo: " + jmohep[2 * i] + " " + jmohep[2 * i + 1] +
+             " jda: " + jdahep[2 * i] + " " + jdahep[2 * i + 1]);
+             System.out.println("p: " + phep[5 * i] + " " + phep[5 * i + 1] + " " +
+             phep[5 * i + 2] + " " + phep[5 * i + 3] + " " + phep[5 * i + 4]);
+             System.out.println("v: " + vhep[4 * i] + " " + vhep[4 * i + 1] + " " +
+             vhep[4 * i + 2] + " " + vhep[4 * i + 3]);
              */
             //                                }
             i++;
@@ -433,20 +555,23 @@
         sw.writeRecord(ev);
     }
 
-    static private void writeLHEEvent(StreamTokenizer tok, double[] beam,  int nevhep) throws IOException {
-         Random generator = new Random();
+    static private void writeLHEEvent(StreamTokenizer tok, double[] beam, int nevhep) throws IOException {
+        Random generator = new Random();
         getToNextEvent(tok);
         List<Double> nums = getNumbersInLine(tok);
 
-        if (nums.size() != 6)
+        if (nums.size() != 6) {
             throw new RuntimeException("Unexpected entry for number of particles");
+        }
         int nhep = nums.get(0).intValue();
 //        System.out.println("Number of particles for event " + nevhep + ": " + nhep);
 
         double decLen = 0;
         double maxWght = 0;
-          
-        if (expDecay) maxWght = findMaxWeight();
+
+        if (expDecay) {
+            maxWght = findMaxWeight();
+        }
         int isthep[] = new int[nhep];
         int idhep[] = new int[nhep];
         int jmohep[] = new int[2 * nhep];
@@ -456,12 +581,15 @@
         int i = 0;
         int idhepTmp = 0;
         double[] ApMom = {0, 0, 0};
-        double ApMass=0;
-        double ApEnergy=0;
+        double ApMass = 0;
+        double ApEnergy = 0;
+        boolean acceptEvent = false;
+         int nElePass = 0;
         for (int npart = 0; npart < nhep; npart++) {
             List<Double> vals = getNumbersInLine(tok);
-            if (vals.size() != 13)
+            if (vals.size() != 13) {
                 throw new RuntimeException("Unexpected entry for a particle");
+            }
             if (vals.get(1).intValue() != 9) {//ignore the vertex for now
 //                int ip = vals.get(0).intValue();
 //                if (ip != i + 1) {
@@ -469,38 +597,56 @@
 //                }
                 isthep[i] = vals.get(1).intValue();
 
-                if (vals.get(1).intValue() == -1)
+                if (vals.get(1).intValue() == -1) {
                     isthep[i] = 3;
+                }
                 idhepTmp = vals.get(0).intValue();
 //                System.out.println(idhepTmp);
                 idhep[i] = idhepTmp;
-                if (idhep[i] == 611)
+                if (idhep[i] == 611) {
                     idhep[i] = 11;
-                if (idhep[i] == -611)
+                }
+                if (idhep[i] == -611) {
                     idhep[i] = -11;
+                }
                 jmohep[2 * i] = vals.get(2).intValue();
                 jmohep[2 * i + 1] = vals.get(3).intValue();
                 jdahep[2 * i] = vals.get(4).intValue();
                 jdahep[2 * i + 1] = vals.get(5).intValue();
-                for (int j = 0; j < 5; j++)
+                for (int j = 0; j < 5; j++) {
                     phep[5 * i + j] = vals.get(j + 6);
+                }
+                Hep3Vector pRot = rotateToDetector(phep[5 * i + 0], phep[5 * i + 1], phep[5 * i + 2]);
+                phep[5 * i + 0] = pRot.x();
+                phep[5 * i + 1] = pRot.y();
+                phep[5 * i + 2] = pRot.z();
                 if (Math.abs(idhepTmp) == 622) {
-                    ApMom[0] = phep[5 * i + 2];
-                    ApMom[1] = phep[5 * i + 1];
-                    ApMom[2] = phep[5 * i + 0];
-                    ApMass=vals.get(11);
-                    ApEnergy=vals.get(10);
+                    /*                  
+                     ApMom[0] = phep[5 * i + 2];
+                     ApMom[1] = phep[5 * i + 1];
+                     ApMom[2] = phep[5 * i + 0];
+                     */
+                    ApMom[0] = pRot.x();
+                    ApMom[1] = pRot.y();
+                    ApMom[2] = pRot.z();
+                    ApMass = vals.get(11);
+                    ApEnergy = vals.get(10);
                 }
-                for (int j = 0; j < 4; j++)
+                for (int j = 0; j < 4; j++) {
                     vhep[4 * i + j] = beam[j];
+                }
 
-                if (!trident && (Math.abs(idhepTmp) == 611 || Math.abs(idhepTmp) == 622)){
-                    double gamma = ApEnergy/ApMass;
-                    if (expDecay) decLen = getDecayLength(maxWght, gamma);
-                     if (flatDecay) decLen = generator.nextDouble() * maxLen;
-                    if (!expDecay && !flatDecay)
+                if (!trident && (Math.abs(idhepTmp) == 611 || Math.abs(idhepTmp) == 622)) {
+                    double gamma = ApEnergy / ApMass;
+                    if (expDecay) {
+                        decLen = getDecayLength(maxWght, gamma);
+                    }
+                    if (flatDecay) {
+                        decLen = generator.nextDouble() * maxLen;
+                    }
+                    if (!expDecay && !flatDecay) {
                         vhep[4 * i + 0] = beam[0] + decLen;
-                    else {                    
+                    } else {
                         double totApMom = Math.sqrt(ApMom[0] * ApMom[0] + ApMom[1] * ApMom[1] + ApMom[2] * ApMom[2]);
                         System.out.println("Decay at : " + decLen);
                         vhep[4 * i + 0] = beam[0] + decLen * ApMom[0] / totApMom;
@@ -508,15 +654,23 @@
                         vhep[4 * i + 2] = beam[2] + decLen * ApMom[2] / totApMom;
                     }
                 }
+
+                if (Math.abs(idhepTmp) == 611 || Math.abs(idhepTmp) == 11) {
+                    if (inAcceptance(pRot)) {
+//                        acceptEvent = true;
+                        nElePass++;
+                    }
+                }
+
                 // swap x and z axes...
                 /*Don't do this anymore!  We do stuff in JLAB frame now...
-                double px = phep[5 * i + 0];
-                double pz = phep[5 * i + 2];
-                phep[5 * i + 0] = pz;
-                phep[5 * i + 2] = px;
+                 double px = phep[5 * i + 0];
+                 double pz = phep[5 * i + 2];
+                 phep[5 * i + 0] = pz;
+                 phep[5 * i + 2] = px;
 
                  */
-            //                                if (i == 0 || i == nhep - 1) {
+                //                                if (i == 0 || i == nhep - 1) {
 //                    System.out.println(i + " st: " + isthep[i] + " id: " + idhep[i] +
 //                            " jmo: " + jmohep[2 * i] + " " + jmohep[2 * i + 1] +
 //                            " jda: " + jdahep[2 * i] + " " + jdahep[2 * i + 1]);
@@ -525,10 +679,46 @@
 //                    System.out.println("v: " + vhep[4 * i] + " " + vhep[4 * i + 1] + " " +
 //                            vhep[4 * i + 2] + " " + vhep[4 * i + 3]);
 //                                }
-                    i++;
+                i++;
             }
         }
         StdhepEvent ev = new StdhepEvent(nevhep, nhep, isthep, idhep, jmohep, jdahep, phep, vhep);
-        sw.writeRecord(ev);
+        if(nElePass>=_nEleRequired)
+            acceptEvent=true;
+        if (acceptEvent) {
+            sw.writeRecord(ev);
+        } else{
+//            System.out.println("LHE Event failed acceptance cuts");
+        }
+    }
+
+    public static Hep3Vector rotateToDetector(double x, double y, double z) {
+        Hep3Vector vecOrig =
+                new BasicHep3Vector(x, y, z);
+        return rot.rotated(vecOrig);
+    }
+
+    public static String convertDecimal(String num) {
+        if (num.contains(".")) {
+            num = num.replace(".", "pt");
+        }
+        return num;
+    }
+
+    public static boolean inAcceptance(Hep3Vector ph) {
+        boolean ok = false;
+        double[] p = {ph.x(), ph.y(), ph.z()};
+        double ptot = Math.sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
+         
+        double sinThx = p[0] / Math.sqrt(p[0] * p[0] + p[2] * p[2]);
+        double sinThy = p[1] / Math.sqrt(p[1] * p[1] + p[2] * p[2]);
+        // for now, just use thetay (non-bend direction)
+//         System.out.println("px = "+p[0]+"; py = "+p[1]+"; pz = "+p[2]);
+//        System.out.println(sinThy+" "+sinThx+" "+ptot);
+                    
+        if (Math.abs(sinThy) > 0.012&&ptot>0.1) {
+            ok = true;
+        }
+        return ok;
     }
 }

hps-java/src/main/resources/org/lcsim/hps/recon/tracking/strategies
HPS-Test-Lyr50.xml added at 1.1
diff -N HPS-Test-Lyr50.xml
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ HPS-Test-Lyr50.xml	29 Aug 2012 15:40:46 -0000	1.1
@@ -0,0 +1,30 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<StrategyList xmlns:xs="http://www.w3.org/2001/XMLSchema-instance" xs:noNamespaceSchemaLocation="http://lcsim.org/recon/tracking/seedtracker/strategybuilder/strategies.xsd">
+    <TargetDetector>HPS-Test-Lyr50</TargetDetector>
+    <Strategy name="HelicalTrackHit Strategy">
+        
+        <!--Cutoffs-->
+        
+        <MinPT>0.100</MinPT>
+        <MinHits>5</MinHits>
+        <MinConfirm>1</MinConfirm>
+        
+        <MaxDCA>4.0</MaxDCA>
+        <MaxZ0>4.0</MaxZ0>
+
+        <MaxChisq>25.0</MaxChisq>
+        <BadHitChisq>10.0</BadHitChisq>
+
+        <!--Layers-->
+
+        <Layers>                 
+            <Layer type="Seed" layer_number="5" detector_name="Tracker" be_flag="BARREL" />
+            <Layer type="Seed" layer_number="3" detector_name="Tracker" be_flag="BARREL" />
+            <Layer type="Seed" layer_number="1" detector_name="Tracker" be_flag="BARREL" />
+            <Layer type="Confirm" layer_number="7" detector_name="Tracker" be_flag="BARREL" />
+            <Layer type="Extend" layer_number="9" detector_name="Tracker" be_flag="BARREL" />
+             <Layer type="Extend" layer_number="11" detector_name="Tracker" be_flag="BARREL" />            
+        </Layers>
+    </Strategy>
+</StrategyList>
+
CVSspam 0.2.12


Use REPLY-ALL to reply to list

To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1