hps-java/src/test/java/org/lcsim/hps/recon/tracking
diff -u -r1.5 -r1.6
--- TestRunTrackReconTest.java 13 Jul 2013 01:44:56 -0000 1.5
+++ TestRunTrackReconTest.java 3 Aug 2013 00:57:03 -0000 1.6
@@ -15,6 +15,7 @@
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
+import static junit.framework.Assert.assertTrue;
import junit.framework.TestCase;
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
@@ -32,7 +33,7 @@
* @author phansson
*/
public class TestRunTrackReconTest extends TestCase {
-
+
static final String testURLBase = "http://www.slac.stanford.edu/~phansson/files/hps_java_test/HPSTestRunv3/hps-java-1.7-SNAPSHOT-050113";
static final String testFileName = "egs_5.5gev_0.016x0_500mb_recoil_recon_1_hpsTestRunTrackingTest.slcio";
static final String testURLBaseCmp = "http://www.slac.stanford.edu/~phansson/files/hps_java_test/HPSTestRunv3/hps-java-1.7-SNAPSHOT-050113";
@@ -41,15 +42,15 @@
static final boolean saveForReference = false;
static final boolean cmpHistograms = true;
private final int nEvents = 5000;
-
+
public void testTrackRecon() throws Exception {
-
+
File lcioInputFile = null;
-
+
URL testURL = new URL(testURLBase + "/" + testFileName);
FileCache cache = new FileCache();
- lcioInputFile = cache.getCachedFile(testURL);
-
+ lcioInputFile = cache.getCachedFile(testURL);
+
//Process and write out the file
LCSimLoop loop = new LCSimLoop();
loop.setLCIORecordSource(lcioInputFile);
@@ -59,7 +60,7 @@
loop.add(new LCIODriver(outputFile));
loop.loop(nEvents, null);
loop.dispose();
-
+
//Read LCIO back and test!
LCSimLoop readLoop = new LCSimLoop();
readLoop.add(new TrackReconTestDriver());
@@ -67,32 +68,26 @@
readLoop.loop(nEvents, null);
readLoop.dispose();
}
-
-
+
/*
* static nested class that runs the drivers for standard HPS track recon
*/
-
-
-
static class TrackReconTestDriver extends Driver {
private AIDA aida = AIDA.defaultInstance();
private IAnalysisFactory af = aida.analysisFactory();
- private final IHistogram1D hntracks = aida.histogram1D("hntracks",10,0,10);
- private final IHistogram1D hnstereohits = aida.histogram1D("hnstereohits",20,0,20);
- private final IHistogram1D hhitpositionx = aida.histogram1D("hhitpositionx",50,0,800);
- private final IHistogram1D hhitpositiony = aida.histogram1D("hhitpositiony",50,-100,100);
- private final IHistogram1D hhitpositionz = aida.histogram1D("hhitpositionz",50,-100,100);
- private final List<String> histograms = Arrays.asList("hntracks","hnstereohits","hhitpositionx","hhitpositiony","hhitpositionz");
+ private final IHistogram1D hntracks = aida.histogram1D("hntracks", 10, 0, 10);
+ private final IHistogram1D hnstereohits = aida.histogram1D("hnstereohits", 20, 0, 20);
+ private final IHistogram1D hhitpositionx = aida.histogram1D("hhitpositionx", 50, 0, 800);
+ private final IHistogram1D hhitpositiony = aida.histogram1D("hhitpositiony", 50, -100, 100);
+ private final IHistogram1D hhitpositionz = aida.histogram1D("hhitpositionz", 50, -100, 100);
+ private final List<String> histograms = Arrays.asList("hntracks", "hnstereohits", "hhitpositionx", "hhitpositiony", "hhitpositionz");
static final double alpha = 0.05; // Type-I error rate
private int ntracks = 0;
private int nevents = 0;
-
//Test thresholds
private final double ftracks_thr = 0.1;
-
-
+
@Override
protected void process(EventHeader event) {
super.process(event);
@@ -100,25 +95,25 @@
//System.out.printf("%s: found %d tracks\n",this.getName(),tracks.size());
hntracks.fill(tracks.size());
ntracks += tracks.size();
- for(Track track : tracks) {
+ for (Track track : tracks) {
List<TrackerHit> hitsOnTrack = track.getTrackerHits();
hnstereohits.fill(hitsOnTrack.size());
- for(TrackerHit hit : hitsOnTrack) {
+ for (TrackerHit hit : hitsOnTrack) {
double pos[] = hit.getPosition();
this.hhitpositionx.fill(pos[0]);
this.hhitpositiony.fill(pos[1]);
this.hhitpositionz.fill(pos[2]);
}
}
-
+
++nevents;
}
@Override
protected void endOfData() {
super.endOfData();
-
- if(saveForReference) {
+
+ if (saveForReference) {
File outputFile = new TestOutputFile(testFileName.replaceAll(".slcio", "") + "_hpsTestRunTrackingTest.aida");
try {
aida.saveAs(outputFile);
@@ -126,25 +121,25 @@
Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex);
}
}
-
+
/*
* Basic tests
*/
- assertTrue("Failed to find any tracks",ntracks>0);
- double ftracks = ((double)ntracks)/(double)nevents;
- assertTrue("Failed to reconstruct more than " + this.ftracks_thr + " of tracks/event ("+ftracks+")",ftracks>this.ftracks_thr);
- assertTrue("Failed to find any stereo hits",this.hnstereohits.mean()>0.);
-
-
-
+ assertTrue("Failed to find any tracks", ntracks > 0);
+ double ftracks = ((double) ntracks) / (double) nevents;
+ assertTrue("Failed to reconstruct more than " + this.ftracks_thr + " of tracks/event (" + ftracks + ")", ftracks > this.ftracks_thr);
+ assertTrue("Failed to find any stereo hits", this.hnstereohits.mean() > 0.);
+
+
+
// IPlotter plotter = af.createPlotterFactory().create();
// plotter.createRegions(1, 3, 0);
// plotter.setTitle("Nr of tracks");
// plotter.style().statisticsBoxStyle().setVisible(false);
// plotter.region(0).plot(hntracks);
// plotter.show();
-
- if(cmpHistograms) {
+
+ if (cmpHistograms) {
File aidaCmpInputFile = null;
@@ -154,12 +149,11 @@
FileCache cache = null;
try {
cache = new FileCache();
-aida.saveAs("testResults.aida");
} catch (IOException ex) {
Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex);
}
try {
- aidaCmpInputFile = cache.getCachedFile(cmpURL);
+ aidaCmpInputFile = cache.getCachedFile(cmpURL);
} catch (IOException ex) {
Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex);
}
@@ -171,23 +165,28 @@
} catch (IOException ex) {
Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex);
}
- if(tree_cmp==null) {
+ if (tree_cmp == null) {
throw new RuntimeException("cannot create the ITree for the comparison file located at" + aidaCmpInputFile.getAbsolutePath());
}
- for(String histname : histograms) {
-
- IHistogram1D h_ref = (IHistogram1D)tree_cmp.find(histname);
+ for (String histname : histograms) {
+
+ IHistogram1D h_ref = (IHistogram1D) tree_cmp.find(histname);
IHistogram1D h_test = aida.histogram1D(histname);
- boolean nullHypoIsRejected = CompareHistograms.instance().getTTest(alpha, h_test.mean(), h_ref.mean(), h_test.rms()*h_test.rms(), h_ref.rms()*h_ref.rms(), h_test.allEntries(), h_ref.allEntries());
- double p_value = CompareHistograms.instance().getTTestPValue(h_test.mean(), h_ref.mean(), h_test.rms()*h_test.rms(), h_ref.rms()*h_ref.rms(), h_test.allEntries(), h_ref.allEntries());
+ boolean nullHypoIsRejected = CompareHistograms.instance().getTTest(alpha, h_test.mean(), h_ref.mean(), h_test.rms() * h_test.rms(), h_ref.rms() * h_ref.rms(), h_test.allEntries(), h_ref.allEntries());
+ double p_value = CompareHistograms.instance().getTTestPValue(h_test.mean(), h_ref.mean(), h_test.rms() * h_test.rms(), h_ref.rms() * h_ref.rms(), h_test.allEntries(), h_ref.allEntries());
+ System.out.printf("%s: %s %s T-Test (%.1f%s C.L.) with p-value=%.3f\n", TestRunTrackReconTest.class.getName(), histname, (nullHypoIsRejected ? "FAILED" : "PASSED"), (1 - alpha) * 100, "%", p_value);
+ assertTrue("Failed T-Test (" + (1 - alpha) * 100 + "% C.L. p-value=" + p_value + ") comparing histogram " + histname, !nullHypoIsRejected);
+ double ks_p_value = CompareHistograms.getKolmogorovPValue(h_test, h_ref);
+ boolean ksNullHypoIsRejected = (ks_p_value < alpha);
+ System.out.printf("%s: %s %s Kolmogorov-Smirnov test (%.1f%s C.L.) with p-value=%.3f\n", TestRunTrackReconTest.class.getName(), histname, (ksNullHypoIsRejected ? "FAILED" : "PASSED"), (1 - alpha) * 100, "%", ks_p_value);
+ assertTrue("Failed Kolmogorov-Smirnov test (" + (1 - alpha) * 100 + "% C.L. p-value=" + ks_p_value + ") comparing histogram " + histname, !ksNullHypoIsRejected);
+
//TODO: use a real two-sample Poisson test
// boolean entriesInconsistent = CompareHistograms.instance().getTTest(alpha, h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries());
// double p_value_entries = CompareHistograms.instance().getTTestPValue(h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries());
- System.out.printf("%s: %s %s T-Test (%.1f%s C.L.) with p-value=%.3f\n",TestRunTrackReconTest.class.getName(),histname,(nullHypoIsRejected?"FAILED":"PASSED"),(1-alpha)*100,"%",p_value);
// System.out.printf("%s: %s entries are %s (N=%d Ref=%d for T-Test w/ %.1f%s C.L.)\n",TestRunTrackReconTest.class.getName(),histname,(entriesInconsistent?"INCONSISTENT":"CONSISTENT"),h_test.entries(),h_ref.entries(),(1-alpha)*100,"%");
- assertTrue("Failed T-Test ("+ (1-alpha)*100 + "% C.L. p-value=" + p_value + ") comparing histogram " + histname,!nullHypoIsRejected);
// assertTrue("Failed T-Test ("+ (1-alpha)*100 + "% C.L. p-value=" + p_value_entries + ") on number of entries comparing histogram " + histname,!entriesInconsistent);
-
+
}
@@ -195,13 +194,7 @@
Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex);
}
}
-
- }
-
-
-
+ }
}
-
-
}
hps-java/src/main/java/org/lcsim/hps/util
diff -u -r1.1 -r1.2
--- CompareHistograms.java 28 Nov 2012 06:36:29 -0000 1.1
+++ CompareHistograms.java 3 Aug 2013 00:57:03 -0000 1.2
@@ -1,42 +1,42 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
package org.lcsim.hps.util;
+import hep.aida.IAxis;
+import hep.aida.IHistogram1D;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.apache.commons.math.MathException;
import org.apache.commons.math.stat.descriptive.StatisticalSummaryValues;
import org.apache.commons.math.stat.inference.TTestImpl;
+import org.apache.commons.math3.distribution.KolmogorovSmirnovDistribution;
/**
*
- * Class to do various comparisons of histograms
- * Singleton instance with lazy instantiation
- *
+ * Class to do various comparisons of histograms Singleton instance with lazy
+ * instantiation
+ *
* @author phansson
*/
public class CompareHistograms {
-
+
public static CompareHistograms _instance = null;
TTestImpl tTest;
+
private CompareHistograms() {
- tTest = new TTestImpl();
+ tTest = new TTestImpl();
}
-
+
public static CompareHistograms instance() {
- if(_instance==null) {
+ if (_instance == null) {
_instance = new CompareHistograms();
}
return _instance;
}
-
+
public double getTTestPValue(double m1, double m2, double v1, double v2, int n1, int n2) {
- StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1,n1, 1.,0., 0.);
- StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2,n2, 1.,0., 0.);
+ StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1, n1, 1., 0., 0.);
+ StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2, n2, 1., 0., 0.);
- double p_value = -1;
+ double p_value = -1;
try {
p_value = tTest.tTest(stat1, stat2);
} catch (IllegalArgumentException ex) {
@@ -45,12 +45,12 @@
Logger.getLogger(CompareHistograms.class.getName()).log(Level.SEVERE, null, ex);
}
return p_value;
-
+
}
- public boolean getTTest(double alpha,double m1, double m2, double v1, double v2, int n1, int n2) {
- StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1,n1, 1.,0., 0.);
- StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2,n2, 1.,0., 0.);
+ public boolean getTTest(double alpha, double m1, double m2, double v1, double v2, int n1, int n2) {
+ StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1, n1, 1., 0., 0.);
+ StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2, n2, 1., 0., 0.);
boolean nullHypoIsRejected = false;
try {
@@ -62,8 +62,26 @@
}
return nullHypoIsRejected;
}
-
-
-
+ public static double getKolmogorovPValue(IHistogram1D reference, IHistogram1D test) {
+ double integralMax = 0.0;
+ double refIntegral = reference.binHeight(IAxis.UNDERFLOW_BIN);
+ double testIntegral = test.binHeight(IAxis.UNDERFLOW_BIN);
+ double integralDiff = Math.abs(refIntegral / reference.allEntries() - testIntegral / test.allEntries());
+ if (integralDiff > integralMax) {
+ integralMax = integralDiff;
+ }
+ for (int i = 0; i < reference.axis().bins(); i++) {
+ refIntegral += reference.binHeight(i);
+ testIntegral += test.binHeight(i);
+
+ integralDiff = Math.abs(refIntegral / reference.allEntries() - testIntegral / test.allEntries());
+ if (integralDiff > integralMax) {
+ integralMax = integralDiff;
+ }
+ }
+ int n = (int) Math.ceil(Math.sqrt((reference.allEntries() * test.allEntries()) / (reference.allEntries() + test.allEntries())));
+ KolmogorovSmirnovDistribution dist = new KolmogorovSmirnovDistribution(n);
+ return 1.0 - dist.cdf(integralMax);
+ }
}