hps-java/src/main/java/org/lcsim/HPSUsers/Example
diff -u -r1.1 -r1.2
--- DetailedAnalysisDriver.java 20 Apr 2011 20:07:39 -0000 1.1
+++ DetailedAnalysisDriver.java 19 May 2011 17:21:48 -0000 1.2
@@ -120,7 +120,8 @@
public String outputTextName = "myevents.txt";
FileWriter fw;
PrintWriter pw;
- double[] beamsize = {0.001, 0.02, 0.02};
+ double[] beamsize = {0.001, 0.02, 0.02};
+
public DetailedAnalysisDriver(int layers) {
nlayers[0] = layers;
@@ -280,13 +281,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");
+ if (event.hasItem("AxialTrackHits")) {
+ List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
- int nAxialHitsTotal = axialhits.size();
- int nL1Hits = 0;
- for (HelicalTrackHit hth : axialhits)
- if (hth.Layer() == 1)
- nL1Hits++;
+ int nAxialHitsTotal = axialhits.size();
+ int nL1Hits = 0;
+ for (HelicalTrackHit hth : axialhits)
+ if (hth.Layer() == 1)
+ nL1Hits++;
+ }
Map<String, Integer> occupancyMap = new HashMap<String, Integer>();
for (RawTrackerHit rh : rawHits) {
@@ -340,24 +343,30 @@
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());
+ if (event.hasItem("AxialTrackMCRelations")) {
+ 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");
aida.cloud1D("Matched Tracks per Event").fill(tracklist.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);
-
+ List<Track> axialtracklist = null;
+ RelationalTable trktomcAxial = null;
+ if (event.hasItem("AxialTracks")) {
+ axialtracklist = event.get(Track.class, "AxialTracks");
+ aida.cloud1D("Axial Tracks per Event").fill(axialtracklist.size());
+ trktomcAxial = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ }
int _nchRec = 0;
int _neleRec = 0;
int _nposRec = 0;
@@ -382,11 +391,13 @@
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);
+ if (event.hasItem("AxialTracks")) {
+ 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
@@ -705,13 +716,15 @@
}
//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());
+ if (event.hasItem("AxialTracks")) {
+ 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) {
@@ -801,7 +814,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.getMCParticle();
if (mcp == null)
@@ -817,11 +830,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;
@@ -838,7 +851,7 @@
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;
@@ -858,14 +871,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);
@@ -873,38 +886,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());
@@ -916,15 +929,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());
}
}
}
@@ -970,10 +983,11 @@
// Check cases where we have multiple tracks associated with this MC particle
Set<Track> trklist = trktomc.allTo(mcp);
int ntrk = trklist.size();
-
-
- Set<Track> trklistAxial = trktomcAxial.allTo(mcp);
- int ntrkAxial = trklistAxial.size();
+ int ntrkAxial=0;
+ if (event.hasItem("AxialTracks")) {
+ Set<Track> trklistAxial = trktomcAxial.allTo(mcp);
+ ntrkAxial = trklistAxial.size();
+ }
if (mcp.getPDGID() == 622) {
boolean bothreco = true;