lcsim/src/org/lcsim/contrib/Pelham/Example1
diff -N HistogramAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HistogramAnalysisDriver.java 25 Jul 2008 18:26:20 -0000 1.1
@@ -0,0 +1,199 @@
+/*
+ * AnalysisDriver.java
+ *
+ * Created on February 11, 2008, 11:47 AM
+ *
+ */
+
+package org.lcsim.contrib.Pelham.Example1;
+
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map;
+import java.util.HashMap;
+
+import org.lcsim.contrib.seedtracker.SeedTrack;
+import org.lcsim.contrib.seedtracker.SeedCandidate;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+
+
+/**
+ *
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class HistogramAnalysisDriver extends Driver {
+ private AIDA aida = AIDA.defaultInstance();
+
+
+ /** Creates a new instance of AnalysisDriver */
+ public HistogramAnalysisDriver() {
+
+ }
+
+ /**
+ * Process the current event
+ * @param event EventHeader for this event
+ */
+ public void process(EventHeader event) {
+ List<Track> tracklist = event.getTracks();
+ Map<MCParticle, Track> trkmap = new HashMap<MCParticle, Track>();
+
+ for (Track track : tracklist) {
+ Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ if (track instanceof SeedTrack) {
+ SeedCandidate seed = ((SeedTrack) track).getSeedCandidate();
+ for (HelicalTrackHit hit : seed.getHits()) {
+
+
+
+
+ List<MCParticle> mclist = hit.getMCParticles();
+ for (MCParticle mcp : mclist) {
+ if (mcmap.containsKey(mcp)) {
+ int entries = mcmap.get(mcp) + 1;
+ mcmap.put(mcp, entries);
+ } else
+ mcmap.put(mcp, 1);
+ }
+ }
+ MCParticle mcmatch = null;
+ int nmatch = 0;
+ for (Map.Entry<MCParticle, Integer> me : mcmap.entrySet()) {
+ if (me.getValue() > nmatch) {
+ nmatch = me.getValue();
+ mcmatch = me.getKey();
+ }
+ }
+ if (trkmap.containsKey(mcmatch)) System.out.println("more than one track associated with an MCParticle!!");
+ else trkmap.put(mcmatch, track);
+ }
+ }
+ //find mc particle is creating hits, % hits by majority mc particle
+ //sort hits with most mc particles
+ List<MCParticle> mclist = event.getMCParticles();
+ for (MCParticle mcp : mclist) {
+ if (mcp.getCharge() == 0) continue;
+ if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) continue;
+
+ if (trkmap.containsKey(mcp)) {
+ Track trk = trkmap.get(mcp);
+ if (trk instanceof SeedTrack) {
+ SeedTrack strk = (SeedTrack) trk;
+ HelicalTrackFit fit = strk.getSeedCandidate().getHelix();
+
+ //Class to draw histrogram
+ HelixParamHistograms draw = new HelixParamHistograms(fit,mcp,event);
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, 5);
+
+ //---Helix Parameter Histrograms------------Helix Parameter Histrograms---\\
+
+ if(mcp.getCharge()>0)
+ {
+ draw.DrawResidualPositive();
+ draw.DrawPullPositive();
+ }
+ else
+ {
+ draw.DrawResidualNegative();
+ draw.DrawPullNegative();
+ }
+ draw.DrawResidual();
+ draw.DrawPull();
+
+
+ List<HelicalTrackHit> hit = strk.getSeedCandidate().getHits();
+ TrackHitHistograms draw2 = new TrackHitHistograms(fit, mcp,event, hit);
+
+ //---MCParticle vs HelicalTrackFit----------------------------MCParticle vs HelicalTrackFit---
+ draw2.drawMCParticlevHelicalTrackHit();
+
+ //---SimTracker vs HelicalTrackFit---------------------------SimTracker vs HelicalTrackFit---
+ draw2.drawSimTrackervHelicalTrack();
+
+ Map<MCParticle,Integer> occurrences = new HashMap<MCParticle,Integer>();
+ double purity=0;
+ double mostoccured;
+ int index;
+ for (HelicalTrackHit h : strk.getSeedCandidate().getHits()){
+ List<MCParticle> mcpl = h.getMCParticles();
+ for(MCParticle part : mcpl){
+ if(occurrences.containsKey(part)){
+ index = occurrences.get(part);
+ occurrences.put(part,index+1);
+ }
+ else {
+ occurrences.put(part, 1);
+ }
+ }
+ //(1-purity)*totalhits
+
+ }
+ if(occurrences.size()==1){
+ purity =1;
+ }
+ else{
+ mostoccured = -1;
+ for(MCParticle p : occurrences.keySet()){
+ int occured = occurrences.get(p);
+ if(occured>mostoccured){
+ mostoccured = occured;
+ }
+ }
+ purity = mostoccured / strk.getSeedCandidate().getHits().size();
+ }
+ aida.histogram1D("Purity/Purity", 200, 0, 1.1).fill(purity);
+ aida.histogram2D("Purity/PvsHits", 200, 0, 1.1,200,0,20).fill(purity,strk.getSeedCandidate().getHits().size());
+ if(purity<1){
+ aida.histogram1D("Purity/Momentum", 300, 0, 150).fill(calc.getMCMomentum());
+ aida.histogram1D("Purity/MCPt", 300, 0, 150).fill(calc.getMCTransverseMomentum());
+ aida.histogram1D("Purity/Omega", 300, -.01, .01).fill(calc.getMCOmega());
+ aida.histogram1D("Purity/tanL", 300, -6, 6).fill(calc.getSlopeSZPlane());
+ aida.histogram1D("Purity/z0", 300, -6, 6).fill(calc.getZ0());
+ aida.histogram1D("Purity/1-purity", 200, 0, 10).fill((1-purity)*strk.getSeedCandidate().getHits().size());
+ aida.histogram1D("Purity/Purity to hits: "+strk.getSeedCandidate().getHits().size(), 300, 0, 1).fill(purity);
+ }
+
+
+ }
+ }
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, 5);
+ double wgt=0;
+ double theta = calc.getTheta();
+ double pt = calc.getMCTransverseMomentum();
+ double degreetheta = 180*theta / Math.PI;
+ if(trkmap.containsKey(mcp)){
+ wgt=1;
+ }
+ if(pt>1.1){
+ aida.profile1D("Efficiency/Efficiency vs theta", 90, 0., 180.).fill(degreetheta, wgt);
+ aida.histogram1D("Efficiency/MC angle", 90, 0., 180.).fill(degreetheta);
+ }
+ if (Math.cos(theta) < 0.985) {
+ aida.profile1D("Efficiency/Efficiency vs pT", 100, 0., 50.).fill(pt, wgt);
+ aida.histogram1D("Efficiency/MC pT", 100, 0., 50.).fill(pt);
+ }
+ if (!trkmap.containsKey(mcp)) {
+ List<HelicalTrackHit> hits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
+ if (hits.size() > 6) {
+ System.out.println("Failed to find track. Found "+hits.size()+" hits @ theta = "
+ +180.*Math.acos(mcp.getMomentum().z() / mcp.getMomentum().magnitude())/Math.PI);
+ for (HelicalTrackHit hit : hits) {
+ System.out.println("Hit in "+hit.getLayerIdentifier()+" at x= "+hit.x()+" y= "+hit.y()+" z= "+hit.z());
+ }
+ }
+ }
+ }
+
+ return;
+ }
+ }
+
