lcsim/src/org/lcsim/contrib/uiowa
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;
+ }
+}