diff -N MipQualityDecision.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MipQualityDecision.java	6 Jun 2008 00:07:27 -0000	1.1
@@ -0,0 +1,129 @@

+package org.lcsim.contrib.uiowa;
+
+import java.util.*; 
+
+import org.lcsim.util.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.*;
+import org.lcsim.event.*;
+
+/**
+ * A class to look at whether a calorimeter cluster is
+ * consistent with being a track/MIP. Criteria are based on
+ *
+ * a) number of layers in cluster
+ * b) number of layers with >1 hit
+ * c) chi^2 of helix fit to cluster
+ *
+ * @author [log in to unmask]
+ * @version $Id: MipQualityDecision.java,v 1.1 2008/06/06 00:07:27 mcharles Exp $
+ */
+ 
+
+public class MipQualityDecision implements DecisionMakerSingle<Cluster>
+{
+    public MipQualityDecision() {
+	// No user-settable options for now
+    }
+
+    public boolean valid(Cluster mip) {
+	// How many layers in this MIP?
+	// And how many instances of >1 hit in a layer?
+	int countDuplicates = 0;
+	Set<Integer> layersSeen = new HashSet<Integer>();
+	for (CalorimeterHit hit : mip.getCalorimeterHits()) {
+	    int layer = getLayer(hit);
+	    Integer layerInt = new Integer(layer);
+	    if (!layersSeen.contains(layerInt)) {
+		layersSeen.add(layerInt);
+	    } else {
+		countDuplicates++;
+	    }
+	}
+	int countHits = mip.getCalorimeterHits().size();
+	int countLayers = layersSeen.size();
+	double fractionDuplicates = ((double)(countDuplicates)) / ((double)(countLayers));
+
+	// Now check chisq. Start by getting the hits:
+	List<org.lcsim.fit.helicaltrack.HelicalTrackHit> hitsToFit = new Vector<org.lcsim.fit.helicaltrack.HelicalTrackHit>();
+	for (CalorimeterHit calHit : mip.getCalorimeterHits()) {
+	    double pos[] = calHit.getPosition();
+	    String subdetName = calHit.getSubdetector().getName();
+	    double err_rphi = 0.0;
+	    double err_z = 0.0;
+	    if (subdetName.compareTo("EMBarrel")==0) {
+		// Barrel -- r fixed but phi uncertain 
+		//   rphi error = cell size / sqrt(12)
+		//   z error = cell size / sqrt(12)
+		err_rphi = 3.5 / 3.464;
+		err_z = 3.5 / 3.464;
+	    } else if (subdetName.compareTo("EMEndcap")==0) {
+		// Endcap -- z fixed but r and phi uncertain
+		//   rphi error = cell size / sqrt(12) ??
+		//   z error = silicon width / sqrt(12)
+		err_rphi = 3.5 / 3.464;
+		err_z = 0.32 / 3.464; // 320 micron silicon
+	    } else if (subdetName.compareTo("HADBarrel")==0) {
+		err_rphi = 10.0 / 3.464;
+		err_z = 10.0 / 3.464;
+	    } else if (subdetName.compareTo("HADEndcap")==0) {
+		err_rphi = 10.0 / 3.464;
+		err_z = 3.0 / 3.464; // 1.2 mm gas OR 5.0 mm scint... pick 3mm as a compromise
+	    } else {
+		throw new AssertionError("Subdet not recognized: "+subdetName);
+	    }
+	    org.lcsim.fit.helicaltrack.HelicalTrack3DHit fitHit = new org.lcsim.fit.helicaltrack.HelicalTrack3DHit(pos[0], pos[1], pos[2], err_rphi, err_z);
+	    hitsToFit.add(fitHit);
+	}
+
+	// Now do the fit:
+	org.lcsim.fit.helicaltrack.HelicalTrackFitter fitter = new org.lcsim.fit.helicaltrack.HelicalTrackFitter();
+	org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus fitOK = fitter.fit(hitsToFit);
+
+	// Check fit quality:
+	if (fitOK == org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus.Success) {
+	    org.lcsim.fit.helicaltrack.HelicalTrackFit fitResult = fitter.getFit();
+	    double chisqArray[] = fitResult.chisq();
+	    int ndfArray[] = fitResult.ndf();
+	    if (chisqArray.length != ndfArray.length) {
+		throw new AssertionError("Chisq length = "+chisqArray.length+" but ndf length = "+ndfArray.length);
+	    }
+	    if (ndfArray.length >= 1) {
+		double chisq = chisqArray[0];
+		double ndf = ndfArray[0];
+		if (ndf > 0) {
+		    double prob = 1.0;
+		    if (chisq > 0.0) {
+			prob = org.lcsim.math.chisq.ChisqProb.gammq(ndf,chisq);
+		    }
+		    double quality = prob * (1.0 - fractionDuplicates);
+		    double cut = 0.05; // require 5% prob by default
+		    if (countLayers < 6) {
+			cut = 0.1; // More stringent for short MIPs which are more likely to be combinatorics
+		    } else if (countLayers > 15) {
+			cut = 0.01; // Less stringent for long MIPs which are likely to be sound but may be multiple-scattered
+		    }
+		    boolean qualityPass = (quality >= cut);
+		    return qualityPass;
+		} else {
+		    // Not enough points for chisq check
+		    return false;
+		}
+	    } else {
+		// No fit output
+		return false;
+	    }
+	} else {
+	    // Fit failed
+	    return false;
+	}
+    }
+
+    protected int getLayer(CalorimeterHit hit) {
+        org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+        id.setID(hit.getCellID());
+        int layer = id.getLayer();
+        return layer;
+    }
+}