 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/users/mgraham on MAIN | |||
| DetailedAnalysisDriver.java | +108 | -75 | 1.3 -> 1.4 |
| TrackExtrapolationAnalysis.java | +1 | -1 | 1.4 -> 1.5 |
| FastTrackAnalysisDriver.java | +107 | -57 | 1.7 -> 1.8 |
| MainJASDriver.java | +5 | -2 | 1.1 -> 1.2 |
| +221 | -135 | ||
a few minor usability changes to my local drivers
diff -u -r1.3 -r1.4 --- DetailedAnalysisDriver.java 29 Aug 2012 15:40:46 -0000 1.3 +++ DetailedAnalysisDriver.java 19 Sep 2012 22:56:46 -0000 1.4 @@ -25,6 +25,7 @@
import java.util.Set; import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.tracker.silicon.SiTrackerModule;
import org.lcsim.event.EventHeader; import org.lcsim.event.LCRelation; import org.lcsim.event.MCParticle;
@@ -41,6 +42,7 @@
import org.lcsim.fit.helicaltrack.HelixParamCalculator; import org.lcsim.fit.helicaltrack.HelixUtils; import org.lcsim.fit.helicaltrack.TrackDirection;
+import org.lcsim.geometry.Detector;
import org.lcsim.hps.recon.vertexing.BilliorTrack; import org.lcsim.hps.recon.vertexing.BilliorVertex; import org.lcsim.hps.recon.vertexing.StraightLineTrack;
@@ -53,8 +55,8 @@
import org.lcsim.util.aida.AIDA; /**
- * - * @author mgraham
+ + @author mgraham
*/
public class DetailedAnalysisDriver extends Driver {
@@ -120,10 +122,11 @@
public String outputTextName = "myevents.txt";
FileWriter fw;
PrintWriter pw;
- double[] beamsize = {0.001, 0.02, 0.02};
- public DetailedAnalysisDriver(int layers) {
- nlayers[0] = layers;
+ double[] beamsize = {0.001, 0.02, 0.02};
+ int flipsign = 1;
+ public DetailedAnalysisDriver(int layers) {
+// nlayers[0] = layers;
// Define the efficiency histograms
IHistogramFactory hf = aida.histogramFactory();
@@ -141,14 +144,16 @@
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.);
-*/
+ /*
+ 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);
@@ -182,6 +187,15 @@
}
}
+
+
+ @Override
+ public void detectorChanged(Detector detector) {
+ // Setup default stereo pairings, which should work for even number of
+ // modules.
+ nlayers[0] = detector.getSubdetector("Tracker").getLayering().getLayers().getNumberOfLayers();
+ System.out.println("Setting nlayers to "+nlayers[0]);
+ }
@Override
public void process(
EventHeader event) {
@@ -202,8 +216,13 @@
String debugDir = "debugPlots/";
String occDir = "occupancyPlots/";
// Get the magnetic field
- Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
double bfield = event.getDetector().getFieldMap().getField(IP).y();
+// System.out.println("bfield = " + bfield);
+ if (bfield < 0) {
+// System.out.println("Flipping signs of reconstructed tracks!");
+ flipsign = -1;
+ }
List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, "RawTrackerHitMaker_RawTrackerHits");
List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
@@ -281,7 +300,7 @@
// List<HelicalTrackHit> hthits = event.get(HelicalTrackHit.class, "MatchedHTHits");
// List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "RotatedHelicalTrackHits");
- // List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
+ // List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
// int nAxialHitsTotal = axialhits.size();
int nL1Hits = 0;
@@ -333,7 +352,7 @@
// 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)
@@ -352,10 +371,10 @@
// 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);
@@ -405,9 +424,9 @@
double slopeErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex));
double curveErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.curvatureIndex, HelicalTrackFit.curvatureIndex));
_nchRec++;
- if (track.getCharge() < 0)
+ if (track.getCharge() * flipsign < 0)
_neleRec++;
- if (track.getCharge() > 0)
+ if (track.getCharge() * flipsign > 0)
_nposRec++;
SeedTrack stEle = (SeedTrack) track;
@@ -444,21 +463,21 @@
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);
if (badLayers.size() > 0) {
- System.out.println(badLayers.toString());
- System.out.println("Bad Layer code: " + badLayerEle);
+// System.out.println(badLayers.toString());
+// System.out.println("Bad Layer code: " + badLayerEle);
}
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)
@@ -471,9 +490,9 @@
// Make plots for fake, non-fake, and all tracks
if (purity < 0.5) {
- if (track.getCharge() < 0)
+ if (track.getCharge() * flipsign < 0)
_neleFake++;
- if (track.getCharge() > 0)
+ if (track.getCharge() * flipsign > 0)
_nposFake++;
cthfake.fill(cth, 1.0);
phifake.fill(phi, 1.0);
@@ -482,9 +501,9 @@
fillTrackInfo(trackdir, "fake tracks", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, d0, xoca, yoca, z0);
} else {
- if (track.getCharge() < 0)
+ if (track.getCharge() * flipsign < 0)
_neleTru++;
- if (track.getCharge() > 0)
+ if (track.getCharge() * flipsign > 0)
_nposTru++;
cthfake.fill(cth, 0.0);
phifake.fill(phi, 0.0);
@@ -499,7 +518,13 @@
if (nbadZ == 3)
fillTrackInfo(trackdir, "3 Bad Z-hits", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, d0, xoca, yoca, z0);
-
+// System.out.println("Track info :");
+// System.out.println("\t charge = " + track.getCharge() * flipsign);
+// System.out.println("\t d0 = " + d0);
+// System.out.println("\t z0 = " + z0);
+// System.out.println("\t phi0 = " + phi0);
+// System.out.println("\t slope = " + slope);
+// System.out.println("\t curve = " + curve);
// Now analyze MC Particles on this track
MCParticle mcp = tkanal.getMCParticleNew();
if (mcp != null) {
@@ -536,8 +561,16 @@
aida.histogram1D(trackdir + "phi0 Pull", 200, -8, 8).fill((phi0 - phi0mc) / phi0Err);
aida.histogram1D(trackdir + "slope Pull", 200, -8, 8).fill((slope - slopemc) / slopeErr);
aida.histogram1D(trackdir + "curvature Pull", 200, -8, 8).fill((curve - curvemc) / curveErr);
- -
+// System.out.println("MC Particle info :");
+// System.out.println("\t PDGID = " + mcp.getPDGID());
+// if (mcp.getParents().size() > 0)
+// System.out.println("\t mother = " + mcp.getParents().get(0).getPDGID());
+// System.out.println("\t charge = " + mcp.getCharge());
+// System.out.println("\t d0 = " + d0mc);
+// System.out.println("\t z0 = " + z0mc);
+// System.out.println("\t phi0 = " + phi0mc);
+// System.out.println("\t slope = " + slopemc);
+// System.out.println("\t curve = " + curvemc);
BasicHep3Vector axial = new BasicHep3Vector();
@@ -584,7 +617,7 @@
if (htlayer == 1)
l1DeltaZ.put(track, z - zTr);
- if (purity == 1 && track.getCharge() > 0 && nhits == 10) {
+ if (purity == 1 && track.getCharge() * flipsign > 0 && nhits == 10) {
if (clusterlist.get(0).rawhits().size() == 1 && clusterlist.get(1).rawhits().size() == 1) {
aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta y: Layer " + htlayer).fill(y - yTr);
aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta z: Layer " + htlayer).fill(z - zTr);
@@ -803,7 +836,7 @@
MCParticle posMC = null;
for (Track track : tracklist) {
- TrackAnalysis tkanal = tkanalMap.get(track);
+ TrackAnalysis tkanal = tkanalMap.get(track);
// Calculate purity and make appropriate plots
MCParticle mcp = tkanal.getMCParticleNew();
if (mcp == null)
@@ -819,11 +852,11 @@
double phi = Math.atan2(py, px);
double cth = pz / Math.sqrt(pt * pt + pz * pz);
- SeedTrack stEle = (SeedTrack) track;
+ SeedTrack stEle = (SeedTrack) track;
SeedCandidate seedEle = stEle.getSeedCandidate();
- HelicalTrackFit ht = seedEle.getHelix();
- double doca=ht.dca();
- double[] poca={ht.x0(),ht.y0(),ht.z0()};
+ HelicalTrackFit ht = seedEle.getHelix();
+ double doca = ht.dca();
+ double[] poca = {ht.x0(), ht.y0(), ht.z0()};
if (mcp.getCharge() > 0) {
posID = track;
posMC = mcp;
@@ -840,15 +873,15 @@
String vertex = "Vertexing/";
String selected = "Selection/";
String nhitsTotal = "NumberOfHits/";
- List<BilliorTrack> btlist = new ArrayList<BilliorTrack>();
+ List<BilliorTrack> btlist = new ArrayList<BilliorTrack>();
for (Track track1 : tracklist) {
Track ele = null;
Track pos = null;
- int ch1 = track1.getCharge();
+ 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();
+ int ch2 = track2.getCharge() * flipsign;
if (track1 != track2 && ch1 == -ch2) {
ele = track1;
pos = track2;
@@ -860,14 +893,14 @@
ApCand++;
int nElectron = ele.getTrackerHits().size();
int nPositron = pos.getTrackerHits().size();
- BilliorTrack btEle = btMap.get(ele);
+ BilliorTrack btEle = btMap.get(ele);
BilliorTrack btPos = btMap.get(pos);
btlist.clear();
btlist.add(btEle);
btlist.add(btPos);
- -
+ +
BilliorVertex bvertexUC = new BilliorVertex(bfield);
bvertexUC.doBeamSpotConstraint(false);
@@ -875,38 +908,38 @@
BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getVertexPosition();
BasicMatrix bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
- double invMass=bvertexUC.getInvMass(); - - aida.histogram1D(vertex+"BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0)); - aida.histogram1D(vertex+"BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0)); - aida.histogram1D(vertex+"BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0)); - aida.histogram1D(vertex+"BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq()); - aida.histogram1D(vertex+"BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0))); - aida.histogram1D(vertex+"BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1))); - aida.histogram1D(vertex+"BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2))); - aida.histogram1D(vertex+"BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(bvertexUC.getInvMass());
+ double invMass = bvertexUC.getInvMass(); + + aida.histogram1D(vertex + "BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0)); + aida.histogram1D(vertex + "BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0)); + aida.histogram1D(vertex + "BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0)); + aida.histogram1D(vertex + "BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq()); + aida.histogram1D(vertex + "BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0))); + aida.histogram1D(vertex + "BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1))); + aida.histogram1D(vertex + "BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2))); + aida.histogram1D(vertex + "BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(bvertexUC.getInvMass());
aida.cloud1D(apdir + "e+e- Invariant Mass").fill(invMass);
if (eleMC != null && posMC != null && ele == eleID && pos == posID)
aida.cloud1D(apdir + "Matched A' Invariant Mass").fill(invMass);
-
- BilliorVertex bvertex = new BilliorVertex(bfield);
+ + BilliorVertex bvertex = new BilliorVertex(bfield);
bvertex.doBeamSpotConstraint(true);
bvertex.tryNewFormalism(btlist);
BasicMatrix bvtxPos = (BasicMatrix) bvertex.getVertexPosition();
BasicMatrix bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
- aida.histogram1D(vertex+"BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0)); - aida.histogram1D(vertex+"BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0)); - aida.histogram1D(vertex+"BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0)); - aida.histogram1D(vertex+"BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChiSq()); - aida.histogram1D(vertex+"BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0))); - aida.histogram1D(vertex+"BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1))); - aida.histogram1D(vertex+"BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
+ aida.histogram1D(vertex + "BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0)); + aida.histogram1D(vertex + "BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0)); + aida.histogram1D(vertex + "BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0)); + aida.histogram1D(vertex + "BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChiSq()); + aida.histogram1D(vertex + "BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0))); + aida.histogram1D(vertex + "BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1))); + aida.histogram1D(vertex + "BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
- aida.histogram1D(vertex+"BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
+ aida.histogram1D(vertex + "BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
@@ -918,15 +951,15 @@
BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getVertexPosition();
BasicMatrix bsconvtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
- aida.histogram1D(vertex+"BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0)); - aida.histogram1D(vertex+"BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0)); - aida.histogram1D(vertex+"BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0)); - aida.histogram1D(vertex+"BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChiSq()); - aida.histogram1D(vertex+"BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0))); - aida.histogram1D(vertex+"BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1))); - aida.histogram1D(vertex+"BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
+ aida.histogram1D(vertex + "BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0)); + aida.histogram1D(vertex + "BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0)); + aida.histogram1D(vertex + "BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0)); + aida.histogram1D(vertex + "BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChiSq()); + aida.histogram1D(vertex + "BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0))); + aida.histogram1D(vertex + "BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1))); + aida.histogram1D(vertex + "BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
- aida.histogram1D(vertex+"BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
+ aida.histogram1D(vertex + "BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
}
}
}
@@ -959,7 +992,7 @@
// 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]);
+ boolean isFindable = findable.InnerTrackerIsFindable(mcp, nlayers[0]);
Set<SimTrackerHit> mchitlist = mcHittomcP.allTo(mcp);
Set<HelicalTrackCross> hitlist = hittomc.allTo(mcp);
@@ -988,7 +1021,7 @@
bothreco = false; // if (findable.LayersHit(d) != nlayers[0]) // if (!findable.InnerTrackerIsFindable(d, nlayers[0] - 2))
- if (!findable.InnerTrackerIsFindable(d, nlayers[0]))
+ if (!findable.InnerTrackerIsFindable(d, nlayers[0]))
bothfindable = false;
}
double vtxWgt = 0;
diff -u -r1.4 -r1.5 --- TrackExtrapolationAnalysis.java 9 Sep 2012 02:56:51 -0000 1.4 +++ TrackExtrapolationAnalysis.java 19 Sep 2012 22:56:46 -0000 1.5 @@ -141,7 +141,7 @@
HPSEcalCluster clust = findClosestCluster(posAtEcal, clusters);
if (clust != null) {
- System.out.println("Cluster Position = ("+clust.getPosition()[0]+","+clust.getPosition()[1]+","+clust.getPosition()[2]+")");
+ // System.out.println("Cluster Position = ("+clust.getPosition()[0]+","+clust.getPosition()[1]+","+clust.getPosition()[2]+")");
double minFringe=HPSTrack.DIPOLE_EDGE-10;
double maxFringe=HPSTrack.DIPOLE_EDGE+50;
// Hep3Vector posAtEcalHPS = hpstrk.getPositionAtZ(zCluster, minFringe,maxFringe, 0.1);
diff -u -r1.7 -r1.8 --- FastTrackAnalysisDriver.java 5 Sep 2012 19:40:05 -0000 1.7 +++ FastTrackAnalysisDriver.java 19 Sep 2012 22:56:46 -0000 1.8 @@ -1,5 +1,6 @@
package org.lcsim.hps.users.mgraham;
+import hep.aida.IAnalysisFactory;
import hep.physics.matrix.BasicMatrix; import hep.physics.vec.BasicHep3Vector; import hep.physics.vec.Hep3Vector;
@@ -40,10 +41,11 @@
import org.lcsim.recon.tracking.seedtracker.SeedCandidate; import org.lcsim.recon.tracking.seedtracker.SeedTrack; import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
/**
- * - * @author mgraham
+ + @author mgraham
*/
public class FastTrackAnalysisDriver extends Driver {
@@ -63,12 +65,14 @@
// 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;
+ boolean makePlots = true; + private AIDA aida = AIDA.defaultInstance(); + private IAnalysisFactory af = aida.analysisFactory();
public FastTrackAnalysisDriver(int trackerLayers, int mintrkLayers, String config) {
nlayers[0] = trackerLayers;
_minLayers = mintrkLayers;
_config = config;
-
}
public FastTrackAnalysisDriver() {
@@ -82,6 +86,7 @@
_minLayers = mintrkLayers;
}
+
public void setBeamSigmaX(double sigma){
beamsize[1]=sigma; //the beamsize[] array is in tracking frame
}
@@ -109,20 +114,20 @@
Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
double bfield = event.getDetector().getFieldMap().getField(IP).y();
- if(bfield<0){
- flipSign=-1;
+ if (bfield < 0) {
+ flipSign = -1;
}
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, "RotatedHelicalTrackHits");
+ List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "RotatedHelicalTrackHits");
// List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
// List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
int nAxialHitsTotal = stripHits.size();
int nL1Hits = 0;
for (SiTrackerHitStrip1D str : stripHits) {
- if (str.getRawHits().get(0).getLayerNumber()== 1) {
+ if (str.getRawHits().get(0).getLayerNumber() == 1) {
nL1Hits++;
}
}
@@ -130,23 +135,24 @@
// 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, "RotatedMCRelations");
+ 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);
@@ -468,7 +474,7 @@
TrackAnalysis tkanalPos = tkanalMap.get(pos); // Integer nElectron = tkanalEle.getNHitsNew(); // Integer nPositron = tkanalPos.getNHitsNew();
- List<Integer> layersEle = tkanalEle.getTrackLayerList();
+ List<Integer> layersEle = tkanalEle.getTrackLayerList();
List<Integer> layersPos = tkanalPos.getTrackLayerList();
Integer nElectron = encodeLayers(layersEle);
Integer nPositron = encodeLayers(layersPos);
@@ -615,6 +621,52 @@
// print out some event information
pw.format("%d %d %d %d", _neleRec, _nposRec, nAxialHitsTotal, nL1Hits);
pw.println();
+ if (makePlots) {
+ aida.histogram1D("Event Number", 50, 0, 100000).fill(nevt);
+ aida.histogram1D("pxE", 50, 0, 2.5).fill(pxE);
+ aida.histogram1D("pyE", 50, -0.1, 0.1).fill(pyE);
+ aida.histogram1D("pzE", 50, -0.1, 0.1).fill(pzE);
+ aida.histogram1D("d0E", 50, -1.0, 1.0).fill(d0E);
+ aida.histogram1D("z0E", 50, -1.0, 1.0).fill(z0E);
+ aida.histogram1D("slopeE", 50, -0.1, 0.1).fill(slopeE);
+ aida.histogram1D("sin(phi0E)", 50, -0.2, 0.2).fill(Math.sin(phi0E));
+ aida.histogram1D("RE", 50, 0, 20000).fill(RE);
+ aida.histogram1D("chisqE", 50, 0, 25).fill(chisqE);
+ aida.histogram1D("l1minE", 50, -10, 10).fill(l1minE);
+ aida.histogram1D("l1dzE", 50, -0.1, 0.1).fill(l1dzE);
+ aida.histogram1D("l123KinkE", 50, -0.1, 0.1).fill(l123KinkE);
+ aida.histogram1D("eleFromAp", 2, 0, 2).fill(eleFromAp);
+
+
+ aida.histogram1D("pxP", 50, 0, 2.5).fill(pxP);
+ aida.histogram1D("pyP", 50, -0.1, 0.1).fill(pyP);
+ aida.histogram1D("pzP", 50, -0.1, 0.1).fill(pzP);
+ aida.histogram1D("d0P", 50, -1.0, 1.0).fill(d0P);
+ aida.histogram1D("z0P", 50, -1.0, 1.0).fill(z0P);
+ aida.histogram1D("slopeP", 50, -0.1, 0.1).fill(slopeP);
+ aida.histogram1D("sin(phi0P)", 50, -0.2, 0.2).fill(Math.sin(phi0P));
+ aida.histogram1D("RP", 50, 0, 20000).fill(RP);
+ aida.histogram1D("chisqP", 50, 0, 25).fill(chisqP);
+ aida.histogram1D("l1minP", 50, -10, 10).fill(l1minP);
+ aida.histogram1D("l1dzP", 50, -10, 10).fill(l1dzP);
+ aida.histogram1D("l123KinkP", 50, -0.1, 0.1).fill(l123KinkP);
+ aida.histogram1D("posFromAp", 2, 0, 2).fill(posFromAp);
+
+ aida.histogram1D("vtxpxE", 50, 0, 2.5).fill(vtxpxE);
+ aida.histogram1D("vtxpyE", 50, -0.1, 0.1).fill(vtxpyE);
+ aida.histogram1D("vtxpzE", 50, -0.1, 0.1).fill(vtxpzE);
+ aida.histogram1D("vtxpxP", 50, 0, 2.5).fill(vtxpxP);
+ aida.histogram1D("vtxpyP", 50, -0.1, 0.1).fill(vtxpyP);
+ aida.histogram1D("vtxpzP", 50, -0.1, 0.1).fill(vtxpzP);
+
+ aida.histogram1D("vtxX", 50, -20, 20).fill(vtx.e(0, 0));
+ aida.histogram1D("vtxY", 50, -0.5, 0.5).fill(vtx.e(1, 0));
+ aida.histogram1D("vtxZ", 50, -0.5, 0.5).fill(vtx.e(2, 0));
+
+ aida.histogram1D("pxmcE", 50, 0, 2.5).fill(pmcEle[0]);
+ aida.histogram1D("pymcE", 50, -0.1, 0.1).fill(pmcEle[1]);
+ aida.histogram1D("pzmcE", 50, -0.1, 0.1).fill(pmcEle[2]);
+ }
}
}
@@ -650,44 +702,44 @@
double pyMCA = -99;
double pzMCA = -99;
double mMCA = -99;
- Hep3Vector decayMCA=new BasicHep3Vector();
+ Hep3Vector decayMCA = new BasicHep3Vector();
int findableE = 0;
int findableP = 0;
int foundE = 0;
int foundP = 0;
- if(mcEle!=null){
- pxMCE = mcEle.getPX();
- 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)) {
- findableP = 1;
- }
- Set<Track> trklistP = trktomc.allTo(mcPos);
- foundP = trklistP.size();//can be greater than 1 if more than 1 track shares hits
- }else{
+ if (mcEle != null) {
+ pxMCE = mcEle.getPX();
+ 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)) {
+ findableP = 1;
+ }
+ Set<Track> trklistP = trktomc.allTo(mcPos);
+ foundP = trklistP.size();//can be greater than 1 if more than 1 track shares hits
+ } 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 !!!!!!!!!");
- }
+ }
+ 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(
@@ -911,12 +963,12 @@
Hep3Vector strvec = VecOp.add(strorigin, struvec);
int strlayer = str.layer();
//for some reason (str!=cl) not working anymore...so make sure the distance between hits isn't 0
- if (layer == strlayer && VecOp.sub(clvec, strvec).magnitude() < Math.abs(mindist)&&VecOp.sub(clvec, strvec).magnitude()!=0) {
+ if (layer == strlayer && VecOp.sub(clvec, strvec).magnitude() < Math.abs(mindist) && VecOp.sub(clvec, strvec).magnitude() != 0) {
mindist = VecOp.sub(clvec, strvec).magnitude();
if (Math.abs(clvec.z()) > Math.abs(strvec.z())) {
mindist = -mindist;
}
- nearest = str;
+ nearest = str;
}
}
@@ -925,9 +977,7 @@
return mindist;
}
-
-
- private Integer encodeLayers(List<Integer> trackLayers) {
+ private Integer encodeLayers(List<Integer> trackLayers) {
Integer hitsEncoded = 0;
for (Integer layer : trackLayers) {
hitsEncoded += (int) Math.pow(2, layer - 1);
diff -u -r1.1 -r1.2 --- MainJASDriver.java 29 Aug 2012 15:40:46 -0000 1.1 +++ MainJASDriver.java 19 Sep 2012 22:56:46 -0000 1.2 @@ -21,9 +21,12 @@
TrackerReconDriver trd=new TrackerReconDriver();
trd.setStripMaxSeparation(20.0);
trd.setStripTolerance(1.0);
- trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-Test-4pt0.xml");
+ trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-Full.xml");
+// trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-Test-Lyr50.xml");
add(trd);
- add(new DetailedAnalysisDriver(10));
+ add(new DetailedAnalysisDriver(12));
+// add(new FastTrackAnalysisDriver());
+
}
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