Commit in hps-java/src/main on MAIN | |||
java/org/lcsim/hps/users/mgraham/MainJASDriver.java | +31 | added 1.1 | |
/FastTrackAnalysisDriver.java | +178 | -89 | 1.2 -> 1.3 |
/DetailedAnalysisDriver.java | +68 | -64 | 1.2 -> 1.3 |
java/org/lcsim/hps/recon/tracking/FindableTrack.java | +1 | -1 | 1.1 -> 1.2 |
resources/org/lcsim/hps/steering/HPSTrackingBunchedMC.lcsim | +58 | added 1.1 | |
java/org/lcsim/hps/recon/vertexing/BilliorVertex.java | +17 | -14 | 1.2 -> 1.3 |
java/org/lcsim/hps/util/ConvertToStdhep.java | +310 | -120 | 1.4 -> 1.5 |
resources/org/lcsim/hps/recon/tracking/strategies/HPS-Test-Lyr50.xml | +30 | added 1.1 | |
+693 | -288 |
Some new stuff/changes for running vertexing simulation again.
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)); + + } + + +}
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; + }
} }
- -
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);
}
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();
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>
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;
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;
} }
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> +
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