hps-java/src/main/java/org/lcsim/hps/examples
diff -u -r1.1 -r1.2
--- StarterAnalysisDriver.java 8 Mar 2012 23:40:51 -0000 1.1
+++ StarterAnalysisDriver.java 8 Feb 2013 02:59:40 -0000 1.2
@@ -1,59 +1,128 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
package org.lcsim.hps.examples;
-import java.io.IOException;
import java.util.List;
-import java.util.logging.Level;
-import java.util.logging.Logger;
import org.lcsim.event.EventHeader;
import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
import org.lcsim.event.RelationalTable;
import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.event.base.BaseRelationalTable;
-import org.lcsim.hps.recon.tracking.TrackAnalysis;
+import org.lcsim.event.base.BaseTrackState;
+import org.lcsim.hps.users.meeg.LCIOTrackAnalysis;
import org.lcsim.util.Driver;
-/**
- *
- * @author mgraham
+/*
+ * Example analysis driver.
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: StarterAnalysisDriver.java,v 1.2 2013/02/08 02:59:40 meeg Exp $
*/
-public class StarterAnalysisDriver extends Driver{
- int nevents=0;
- int naccepted=0;
-
- public void process(
- EventHeader event) {
- nevents++;
+public class StarterAnalysisDriver extends Driver {
+
+ int nevents = 0;
+ int naccepted = 0;
+ private boolean debug = false;
+ private double bfield = 0.5;
+
+ @Override
+ public void process(EventHeader event) {
+ nevents++;
List<Track> tracklist = event.get(Track.class, "MatchedTracks");
- if(tracklist.size()<2)
+ if (tracklist.size() < 2) {
return;
-
+ }
+ System.out.println("Event with "+tracklist.size()+" tracks");
+
+ if (debug) {
+ List<List<TrackerHit>> hitlists = event.get(TrackerHit.class);
+ for (List<TrackerHit> hitlist : hitlists) {
+ System.out.println(event.getMetaData(hitlist).getName());
+ for (TrackerHit hit : hitlist) {
+ System.out.println(hit);
+ }
+ }
+ }
+
RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
-
+ if (debug) {
+ System.out.println("HelicalTrackMCRelations:" + mcrelations.size());
+ }
for (LCRelation relation : mcrelations) {
if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
hittomc.add(relation.getFrom(), relation.getTo());
+ if (debug) {
+ System.out.println("to:" + relation.getTo());
+ System.out.println("from:" + relation.getFrom());
+ }
+ }
+ }
+
+ RelationalTable hittostrip = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> hitrelations = event.get(LCRelation.class, "HelicalTrackHitRelations");
+ if (debug) {
+ System.out.println("HelicalTrackHitRelations:" + hitrelations.size());
+ }
+ for (LCRelation relation : hitrelations) {
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
+ hittostrip.add(relation.getFrom(), relation.getTo());
+ if (debug) {
+ System.out.println("to:" + relation.getTo());
+ System.out.println("from:" + relation.getFrom());
+ }
+ }
+ }
+
+ RelationalTable hittorotated = new BaseRelationalTable(RelationalTable.Mode.ONE_TO_ONE, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> rotaterelations = event.get(LCRelation.class, "RotatedHelicalTrackHitRelations");
+ if (debug) {
+ System.out.println("RotatedHelicalTrackHitRelations:" + rotaterelations.size());
+ }
+ for (LCRelation relation : rotaterelations) {
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
+ hittorotated.add(relation.getFrom(), relation.getTo());
+ if (debug) {
+ System.out.println("to:" + relation.getTo());
+ System.out.println("from:" + relation.getFrom());
+ }
}
}
- int ok=0;
- for (Track track : tracklist) { //remember, these tracks are in the lcsim tracking frame!
- TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
- if(Math.abs(tkanal.getMCParticle().getPDGID())==611)
- ok++;
- //do some stuff to makes sure tracks are great
- //is there an e+e- from the muonium
+
+
+
+ int ok = 0;
+ for (Track track : tracklist) { //remember, these tracks are in the lcsim tracking frame!
+ BaseTrackState ts = (BaseTrackState) track.getTrackStates().get(0);
+ ts.computeMomentum(0.5);
+// BaseTrackState.computeMomentum(track.getTrackStates().get(0), 0.5);
+ LCIOTrackAnalysis tkanal = new LCIOTrackAnalysis(track, hittomc, hittostrip, hittorotated);
+ MCParticle mcp = tkanal.getMCParticle();
+ if (mcp != null) {
+ System.out.println("chisq: " + track.getChi2() + ", pz: " + ts.getMomentum()[0] + ", momentum: " + mcp.getMomentum());
+ if (Math.abs(tkanal.getMCParticle().getPDGID()) == 611) {
+ ok++;
+ }
+ //do some stuff to makes sure tracks are great
+ //is there an e+e- from the muonium
+ }
}
-
- if(ok==2)
+
+ if (ok == 2) {
naccepted++;
-
- }
- public void endOfData() {
-
+ }
+
+ }
+
+ public void setDebug(boolean debug) {
+ this.debug = debug;
+ }
+
+ public void setBfield(double bfield) {
+ this.bfield = bfield;
+ }
+
+ public void endOfData() {
+
System.out.println("# of muonium events= " + naccepted + "; # of total = " + nevents);
}
-}
+}
hps-java/src/main/java/org/lcsim/hps/users/meeg
diff -N LCIOTrackAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LCIOTrackAnalysis.java 8 Feb 2013 02:59:40 -0000 1.1
@@ -0,0 +1,154 @@
+package org.lcsim.hps.users.meeg;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.util.HashMap;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.detector.identifier.Identifier;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.hps.recon.tracking.SvtUtils;
+
+/**
+ *
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: LCIOTrackAnalysis.java,v 1.1 2013/02/08 02:59:40 meeg Exp $
+ */
+public class LCIOTrackAnalysis {
+
+ protected MCParticle _mcp = null;
+ protected double _purity;
+ protected int _nhits;
+ protected int _nbadhits;
+ protected boolean _hasLayerOne;
+ protected Map<Integer, Hep3Vector> _hitLocationPerLayer = new HashMap<Integer, Hep3Vector>();
+ protected int _nhitsNew;
+
+ public LCIOTrackAnalysis(Track trk, RelationalTable hittomc, RelationalTable hittostrip, RelationalTable hittorotated) {
+ // Get the number of hits on the track
+ _nhits = trk.getTrackerHits().size();
+
+ // Create a map containing the number of hits for each MCParticle associated with the track
+ Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ _hasLayerOne = false;
+ // Loop over the hits on the track and make sure we have HelicalTrackHits (which contain the MC particle)
+ for (TrackerHit rotatedHit : trk.getTrackerHits()) {
+ TrackerHit hit = (TrackerHit) hittorotated.from(rotatedHit);
+// System.out.println(hit);
+ // get the set of MCParticles associated with this hit and update the hit count for each MCParticle
+ MCParticle mcp = (MCParticle) hittomc.to(hit);
+// Set<MCParticle> mclist = hittomc.allTo(hit);
+// System.out.println(mclist.size());
+// for (MCParticle mcp : mclist) {
+ if (mcp != null) {
+// System.out.println(mcp.getOrigin());
+ Integer mchits = 0;
+ if (mcmap.containsKey(mcp)) {
+ mchits = mcmap.get(mcp);
+ }
+ mchits++;
+ mcmap.put(mcp, mchits);
+ }
+// }
+
+// BasicHep3Vector axial = new BasicHep3Vector();
+// axial.setV(0, 1, 0);
+
+
+ Set<TrackerHit> hitlist = hittostrip.allTo(hit);
+ for (TrackerHit cl : hitlist) {
+ int layer = SvtUtils.getInstance().getHelper().getValue(new Identifier(cl.getCellID()), "layer");
+ System.out.println(layer);
+ if (layer == 1) {
+ _hasLayerOne = true;
+ }
+
+// _nStripHitsPerLayer[layer - 1] = cl.rawhits().size();
+ _hitLocationPerLayer.put(layer, new BasicHep3Vector(cl.getPosition()));
+ _nhitsNew++;
+// double axdotu = VecOp.dot(cl.u(), axial);
+// boolean isAxial = false;
+// if (axdotu > 0.5) {
+// isAxial = true;
+// _nAxialhits++;
+// } else {
+// _nZhits++;
+// }
+// List<MCParticle> mcPartList = cl.MCParticles();
+// _nMCHitsPerLayer[layer - 1] = mcPartList.size();
+// for (MCParticle mcp : mcPartList) {
+// Integer mchits = 0;
+// if (mcmapAll.containsKey(mcp)) {
+// mchits = mcmapAll.get(mcp);
+// }
+// mchits++;
+// mcmapAll.put(mcp, mchits);
+// if (isAxial) {
+// Integer mchitsAxial = 0;
+// if (mcmapAxial.containsKey(mcp)) {
+// mchitsAxial = mcmapAxial.get(mcp);
+// }
+// mchitsAxial++;
+// mcmapAxial.put(mcp, mchitsAxial);
+// } else {
+// Integer mchitsZ = 0;
+// if (mcmapZ.containsKey(mcp)) {
+// mchitsZ = mcmapZ.get(mcp);
+// }
+// mchitsZ++;
+// mcmapZ.put(mcp, mchitsZ);
+// }
+// }
+ }
+ }
+
+ // Find the MCParticle that has the most hits on the track
+
+ int nbest = 0;
+ MCParticle mcbest = null;
+ for (MCParticle mcp : mcmap.keySet()) {
+ int count = mcmap.get(mcp);
+ if (count > nbest) {
+ nbest = count;
+ mcbest = mcp;
+ }
+ }
+
+ if (nbest > 0) {
+ _mcp = mcbest;
+ }
+ _purity = (double) nbest / (double) _nhits;
+ _nbadhits = _nhits - nbest;
+ }
+
+ public MCParticle getMCParticle() {
+ return _mcp;
+ }
+
+ public int getNHits() {
+ return _nhits;
+ }
+
+ public int getNBadHits() {
+ return _nbadhits;
+ }
+
+ public double getPurity() {
+ return _purity;
+ }
+
+ public int getNHitsNew() {
+ return _nhitsNew;
+ }
+
+ public boolean hasLayerOne() {
+ return _hasLayerOne;
+ }
+
+ public Hep3Vector getClusterPosition(Integer layer) {
+ return _hitLocationPerLayer.get(layer);
+ }
+}