lcsim/src/org/lcsim/contrib/Pelham/Example1
diff -u -r1.2 -r1.3
--- HistogramAnalysisDriver.java 25 Jul 2008 18:28:05 -0000 1.2
+++ HistogramAnalysisDriver.java 7 Aug 2008 21:07:09 -0000 1.3
@@ -8,18 +8,27 @@
package org.lcsim.contrib.Pelham.Example1;
+import hep.aida.IHistogram1D;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Comparator;
import java.util.List;
import java.util.Map;
import java.util.HashMap;
+import java.util.HashSet;
import org.lcsim.contrib.seedtracker.SeedTrack;
import org.lcsim.contrib.seedtracker.SeedCandidate;
+import org.lcsim.contrib.seedtracker.analysis.AnalysisUtils;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
-import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.contrib.seedtracker.analysis.HelixParamCalculator;
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFitter;
+import org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus;
import org.lcsim.util.aida.AIDA;
import org.lcsim.util.Driver;
@@ -30,28 +39,37 @@
* @version 1.0
*/
public class HistogramAnalysisDriver extends Driver {
- private AIDA aida = AIDA.defaultInstance();
-
-
+ private AIDA aida = AIDA.defaultInstance();
+ IHistogram1D thetaeff,pteff,dcaeff,z0eff,phi0eff,tanLeff,omegaeff,eff;
+ private double foundtracks=0,findabletracks=0,effi;
/** Creates a new instance of AnalysisDriver */
public HistogramAnalysisDriver() {
-
- }
+ aida.tree().mkdir("Efficiency");
+ thetaeff = aida.histogramFactory().createHistogram1D("Efficiency/Theta","", 90, 0., 180., "type=efficiency");
+ pteff = aida.histogramFactory().createHistogram1D("Efficiency/Pt","", 100, 0., 50., "type=efficiency");
+ aida.tree().mkdir("Efficiency/HelixParam");
+ dcaeff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/DCA","", 120, -4., 4., "type=efficiency");
+ z0eff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/z0","", 120, -4., 4., "type=efficiency");
+ phi0eff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/phi0","", 120, 0., 360., "type=efficiency");
+ tanLeff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/tanL","", 120, -10., 10., "type=efficiency");
+ omegaeff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/Omega","", 120, -.02, .02, "type=efficiency");
+ eff = aida.histogramFactory().createHistogram1D("Efficiency/efficiency",200,0,3);
+ }
/**
* Process the current event
* @param event EventHeader for this event
*/
- public void process(EventHeader 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()) {
-
+
@@ -76,9 +94,8 @@
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;
@@ -94,7 +111,7 @@
HelixParamCalculator calc = new HelixParamCalculator(mcp, 5);
//---Helix Parameter Histrograms------------Helix Parameter Histrograms---\\
-
+ aida.histogram1D("p", 200, 0, 50).fill(calc.getMCMomentum());
if(mcp.getCharge()>0)
{
draw.DrawResidualPositive();
@@ -110,7 +127,7 @@
List<HelicalTrackHit> hit = strk.getSeedCandidate().getHits();
- TrackHitHistograms draw2 = new TrackHitHistograms(fit, mcp,event, hit);
+ TrackHitHistograms draw2 = new TrackHitHistograms(fit, mcp,event,hit);
//---MCParticle vs HelicalTrackFit----------------------------MCParticle vs HelicalTrackFit---
draw2.drawMCParticlevHelicalTrackHit();
@@ -122,6 +139,7 @@
double purity=0;
double mostoccured;
int index;
+
for (HelicalTrackHit h : strk.getSeedCandidate().getHits()){
List<MCParticle> mcpl = h.getMCParticles();
for(MCParticle part : mcpl){
@@ -133,7 +151,6 @@
occurrences.put(part, 1);
}
}
- //(1-purity)*totalhits
}
if(occurrences.size()==1){
@@ -149,50 +166,178 @@
}
purity = mostoccured / strk.getSeedCandidate().getHits().size();
}
+ int hitsize = strk.getSeedCandidate().getHits().size();
+ String stratname = strk.getStrategy().getName();
+ if(purity<.5)aida.histogram1D("Purity/Purity<.5 only", 200, 0, 1.1).fill(purity);
+ //1D histograms
+ aida.histogram1D("Purity/Strat: "+stratname, 200,0,1.1).fill(purity);
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);
- }
+ if(purity!=1)aida.histogram1D("Purity/Purity<1 only", 200, 0, 1.1).fill(purity);
+ aida.histogram1D("Purity/Purity to hits: "+hitsize, 300, 0, 1.1).fill(purity);
+ //2D histograms
+ aida.histogram1D("Purity/Hits v Purity", 200, 0, 20,"type=efficiency").fill(hitsize,purity);
+ aida.histogram1D("Purity/Theta v Purity", 100, 10, 100, "type=efficiency").fill(180*calc.getTheta()/Math.PI,purity);
+ aida.histogram1D("Purity/Pt v Purity", 200, 0, 70,"type=efficiency").fill(calc.getMCTransverseMomentum(),purity);
-
- }
- }
+
+ }
+ }
HelixParamCalculator calc = new HelixParamCalculator(mcp, 5);
double wgt=0;
double theta = calc.getTheta();
double pt = calc.getMCTransverseMomentum();
+ double dca = calc.getDCA();
+ double z0 = calc.getZ0();
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());
+ Map<MCParticle,List<String>> hitlayermap = new HashMap<MCParticle,List<String>>();
+ Map<MCParticle,List<HelicalTrackHit>> hitordermap = new HashMap<MCParticle,List<HelicalTrackHit>>();
+ //!!
+ List<HelicalTrackHit> listjmk = event.get(HelicalTrackHit.class,"HelicalTrackHits");
+ List<MCParticle> temp = new ArrayList<MCParticle>();
+ AnalysisUtils util = new AnalysisUtils();
+ //!!
+ HelicalTrackFitter fitme = new HelicalTrackFitter();
+
+ Collections.sort(listjmk, new Comparator() {
+
+ public int compare(Object o1, Object o2) {
+ HelicalTrackHit h1 = (HelicalTrackHit) o1;
+ HelicalTrackHit h2 = (HelicalTrackHit) o2;
+
+ double t1, t2;
+ try{
+ t1 = ((RawTrackerHit)h1.getRawHits().get(0)).getSimTrackerHit().get(0).getTime();
+ t2 = ((RawTrackerHit)h2.getRawHits().get(0)).getSimTrackerHit().get(0).getTime();
+ } catch(NullPointerException npe) {
+ t1 = h1.getTime();
+ t2 = h2.getTime();
}
+
+ return Double.compare(t1, t2);
+ }
+ });
+
+ for(HelicalTrackHit hito : listjmk){
+ List<MCParticle> mplist = hito.getMCParticles();
+ for(MCParticle mp : mplist){
+ if(hitlayermap.containsKey(mp)){
+ hitlayermap.get(mp).add(hito.getLayerIdentifier());
+ hitordermap.get(mp).add(hito);
+ }
+ else{
+ ArrayList<String> layers = new ArrayList<String>();
+ layers.add(hito.getLayerIdentifier());
+ hitlayermap.put(mp, layers);
+ ArrayList<HelicalTrackHit> trackhits = new ArrayList<HelicalTrackHit>();
+ trackhits.add(hito);
+ hitordermap.put(mp, trackhits);
+ temp.add(mp);
+ }
+ }
+ }
+ List<MCParticle> testlist = util.FindMCParticleInAllHits(listjmk);
+ if(testlist.containsAll(temp)){System.out.println("true");}
+ HashSet<String> bleh = new HashSet<String>();
+ if(hitlayermap.containsKey(mcp)){
+ for(String e : hitlayermap.get(mcp)){
+ bleh.add(e);
}
}
+
+ if(hitlayermap.containsKey(mcp) && bleh.size()>=7){
+
+ if(trkmap.containsKey(mcp)){
+ wgt=1;
+ eff.fill(1);
+ foundtracks++;
+ findabletracks++;
+ //pass to fitter in sim tracker give any strategy, 'helix fitter" first try helicaltrackfitter
+ }
+ else if(pt>1.1 && Math.abs(dca) < 9 && Math.abs(z0) < 9 && Math.cos(theta) < 0.985 ){
+ findabletracks++;
+ eff.fill(2);
+ if(hitordermap.containsKey(mcp)){
+
+ FitStatus status = fitme.fit(hitordermap.get(mcp));
+ if(status==FitStatus.Success){
+ System.out.println("Fitted Successfully From Analysis Driver");
+ HelicalTrackFit trackfit = fitme.getFit();
+ System.out.println("P: "+mcp.getMomentum().toString());
+ for(HelicalTrackHit h : hitordermap.get(mcp)){
+ System.out.print(h.getLayerIdentifier()+"(time:"+h.getTime()+")");
+ System.out.print("("+h.getPosition()[0]+",");
+ System.out.print(h.getPosition()[1]+",");
+ System.out.print(h.getPosition()[2]+")");
+ System.out.println();
+ }
+ System.out.println("chisq0: "+trackfit.chisq()[0]);
+ System.out.println("chisq1: "+trackfit.chisq()[1]);
+ System.out.println("chisqtot: "+trackfit.chisqtot());
+ System.out.println();
+
+ }
+ else{
+ System.out.println("Failed to fit: "+status.toString());
+ }
+ }
+ aida.histogram1D("Efficiency/NoNHit/MissCount", 100, 0., 2.).fill(1);
+ aida.histogram1D("Efficiency/NoNHit/Size", 100, 0., 15.).fill(hitlayermap.get(mcp).size());
+ aida.histogram1D("Efficiency/NoNHit/p", 100, 0., 10.).fill(calc.getMCMomentum());
+ aida.histogram1D("Efficiency/NoNHit/pt", 100, 0., 14.).fill(calc.getMCTransverseMomentum());
+ //aida.histogram1D("Efficiency/NoNHit/Endpoint", 100, -5000., 5000).fill(mcp.getEndPoint().x()*mcp.getEndPoint().x()+mcp.getEndPoint().y()*mcp.getEndPoint().y());
+ aida.histogram1D("Efficiency/NoNHit/theta", 100, 0., 180).fill(180*calc.getTheta()/Math.PI);
+ aida.histogram1D("Efficiency/NoNHit/theta", 100, 0., 180).fill(180*calc.getTheta()/Math.PI);
+ if(aida.histogram1D("Efficiency/NoNHit/layer hit path/"+hitlayermap.get(mcp).toString(), 100, 0., 360).entries()>1){
+ aida.histogram1D("Efficiency/NoNHit/layer hit path/"+hitlayermap.get(mcp).toString(), 100, 0., 360).setTitle("Interesting");
+ }
+ aida.histogram1D("Efficiency/NoNHit/layer hit path/"+hitlayermap.get(mcp).toString(), 100, 0, 360).fill(180*calc.getTheta()/Math.PI);
+
+ }
+ //particle type
+
+
+ //p vec,endpoint,daughters,charge
+
+ //how many hits
+ //sim final state
+
+ if(pt>1.1 && Math.abs(dca) < 9 && Math.abs(z0) < 9){
+ thetaeff.fill(degreetheta, wgt);
+ tanLeff.fill(calc.getSlopeSZPlane(), wgt);
+ aida.histogram1D("Efficiency/MCType/tanL/"+mcp.getType().getName(),120, -10., 10., "type=efficiency").fill(calc.getSlopeSZPlane(), wgt);
+ }
+ //Pt Cuts: Cos(theta) DCA z0
+ if (Math.cos(theta) < 0.985 && Math.abs(dca) < 9 && Math.abs(z0) < 9) {
+ pteff.fill(pt,wgt);
+ aida.histogram1D("Efficiency/MCType/Pt/"+mcp.getType().getName(), 100, 0., 50., "type=efficiency").fill(pt, wgt);
+ }
+ //Omega Cut: Pt
+ if(pt>1.1){
+ omegaeff.fill(calc.getMCOmega(), wgt);
+ aida.histogram1D("Efficiency/MCType/Omega/"+mcp.getType().getName(), 120, -.02, .02, "type=efficiency").fill(calc.getMCOmega(), wgt);
+ }
+ //Z0 Cuts: Pt DCA Cos(theta)
+ if(pt>1.1 && Math.abs(dca) < 9 && Math.cos(theta) < 0.985){
+ z0eff.fill(calc.getZ0(), wgt);
+ aida.histogram1D("Efficiency/MCType/z0/"+mcp.getType().getName(), 120, -4., 4., "type=efficiency").fill(calc.getZ0(), wgt);
+ }
+ //DCA Cuts: Pt Cos(theta) z0
+ if(pt>1.1 && Math.cos(theta) < 0.985 && Math.abs(z0) < 9){
+ dcaeff.fill(calc.getDCA(), wgt);
+ aida.histogram1D("Efficiency/MCType/DCA/"+mcp.getType().getName(), 120, -4., 4., "type=efficiency").fill(calc.getDCA(), wgt);
+ }
+ //Phi0 Cuts: Pt Cos(theta) z0 DCA
+ if(pt>1.1 && Math.cos(theta) < 0.985 && Math.abs(z0) < 9 && Math.abs(dca) < 9){
+ phi0eff.fill(180*calc.getPhi0()/Math.PI, wgt);
+ aida.histogram1D("Efficiency/MCType/phi0/"+mcp.getType().getName(), 120, 0., 360., "type=efficiency").fill(calc.getDCA(), wgt);
+ }
+ }
}
-
+ effi = (foundtracks/findabletracks)*100;
+ System.out.println("%:"+ effi);
return;
}
+
+
}