diff -N ChargedHadronClusterEnergyCalculator.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ChargedHadronClusterEnergyCalculator.java	19 Nov 2007 21:35:30 -0000	1.1
@@ -0,0 +1,127 @@

+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.conditions.*;
+import org.lcsim.util.*;
+import org.lcsim.event.*;
+import org.lcsim.event.util.*;
+
+/**
+  * ECAL/HCAL calibration designed for charged showers. The hits are divided into two classes: those
+  * which appear to come from MIPs, and all others. The MIP hits are assigned energy based on naive
+  * dE/dx, and the remaining hits are then fed into a ClusterEnergyCalculator as a neutral hadron shower.
+  * Since the number of non-MIP hits could be small or zero, the neutral calibration should not have
+  * a large zero offset.
+  *
+  * The dE/dx energy calculation is currently very crude and hard-coded to the sid01 geometry.
+  * Some improvements needed:
+  *   1) Read in geometry and material information from detector at run-time;
+  *   2) Increase energy loss for tracks that pass through the layer at an angle (currently
+  *      it assumes all tracks are perpendicular to the material planes).
+  *
+  * @version $Id: ChargedHadronClusterEnergyCalculator.java,v 1.1 2007/11/19 21:35:30 mcharles Exp $
+  */
+
+public class ChargedHadronClusterEnergyCalculator extends Driver implements ClusterEnergyCalculator
+{
+    ClusterEnergyCalculator m_neutralCalib = null;
+    String m_mipListName = null;
+    EventHeader m_event = null;
+    public ChargedHadronClusterEnergyCalculator(String mipListName, ClusterEnergyCalculator neutralCalib) {
+	m_neutralCalib = neutralCalib;
+	m_mipListName = mipListName;
+    }
+    public void process(EventHeader event) {
+	m_event = event;
+    }
+    public double getEnergy(Cluster c){
+	// Grab MIP list from the event:
+	List<Cluster> mips = m_event.get(Cluster.class, m_mipListName);
+	Set<CalorimeterHit> mipHits = new HashSet<CalorimeterHit>();
+	for (Cluster mip : mips) {
+	    mipHits.addAll(mip.getCalorimeterHits());
+	}
+	// Create cloned cluster, removing MIP hits
+	double mipEnergySum = 0.0;
+	BasicCluster tmpCopy = new BasicCluster();
+	List<CalorimeterHit> mipHitsInCluster = new Vector<CalorimeterHit>();
+	for (CalorimeterHit hit : c.getCalorimeterHits()) {
+	    if (!mipHits.contains(hit)) {
+		// Not from a MIP segment => treat as neutral hadron
+		tmpCopy.addHit(hit);
+	    } else {
+		// From a MIP segment
+		mipHitsInCluster.add(hit);
+		mipEnergySum += getMIPEnergy(hit);
+	    }
+	}
+	double neutralHadronEnergy = m_neutralCalib.getEnergy(tmpCopy);
+	//System.out.println("DEBUG: Cluster with "+c.getCalorimeterHits().size()+" hits has "+mipHitsInCluster.size()+" mip Hits for an energy of "+neutralHadronEnergy+" [neutral] + "+mipEnergySum+" [mips] = "+(mipEnergySum+neutralHadronEnergy)+". For comparison, naive energy would be "+m_neutralCalib.getEnergy(c));
+	return mipEnergySum+neutralHadronEnergy;
+    }
+    private double getMIPEnergy(Collection<CalorimeterHit> mipHits) {
+	double eSum = 0.0;
+	for (CalorimeterHit hit : mipHits) {
+	    eSum += getMIPEnergy(hit);
+	}
+	return eSum;
+    }
+    // This is ugly for now!
+    private double getMIPEnergy(CalorimeterHit hit) {
+	// First, grab hit info (ID, layer)
+	org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+	id.setID(hit.getCellID());
+	int layer = id.getLayer();
+	// Next, figure out which subdetector we're in.
+	org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+	if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+	String name = subdet.getName();
+	// Set up dE/dx constants:
+	org.lcsim.util.step.DeDx materialDatabase = org.lcsim.util.step.DeDx.instance();
+	double dedx_Silicon = materialDatabase.getDeDx("Silicon");
+	double dedx_Copper = materialDatabase.getDeDx("Copper");
+	double dedx_Kapton = materialDatabase.getDeDx("Kapton");
+	double dedx_Air = materialDatabase.getDeDx("Air");
+	double dedx_Tungsten = materialDatabase.getDeDx("Tungsten");
+	double dedx_Steel235 = materialDatabase.getDeDx("Steel235");
+	double dedx_G10 = materialDatabase.getDeDx("G10");
+	double dedx_RPC = materialDatabase.getDeDx("RPCGasDefault");
+	double dedx_PyrexGlass = materialDatabase.getDeDx("PyrexGlass");
+	// Compute energy loss:
+	double energyLoss = 0.0;
+	if (name.compareTo("EMBarrel")==0 || name.compareTo("EMEndcap")==0) {
+	    if (layer == 0) {
+		energyLoss += 0.032 * dedx_Silicon;
+		energyLoss += 0.005 * dedx_Copper;
+		energyLoss += 0.030 * dedx_Kapton;
+		energyLoss += 0.033 * dedx_Air;
+	    } else if (layer > 0 && layer <= 20) {
+		energyLoss += 0.250 * dedx_Tungsten;
+		energyLoss += 0.032 * dedx_Silicon;
+		energyLoss += 0.005 * dedx_Copper;
+		energyLoss += 0.030 * dedx_Kapton;
+		energyLoss += 0.033 * dedx_Air;
+	    } else if (layer > 20 && layer <= 30) {
+		energyLoss += 0.500 * dedx_Tungsten;
+		energyLoss += 0.032 * dedx_Silicon;
+		energyLoss += 0.005 * dedx_Copper;
+		energyLoss += 0.030 * dedx_Kapton;
+		energyLoss += 0.033 * dedx_Air;
+	    }
+	} else if (name.compareTo("HADBarrel")==0 || name.compareTo("HADEndcap")==0) {
+	    energyLoss += 2.00 * dedx_Steel235;
+	    energyLoss += 0.30 * dedx_G10;
+	    energyLoss += 0.11 * dedx_PyrexGlass;
+	    energyLoss += 0.12 * dedx_RPC;
+	    energyLoss += 0.11 * dedx_PyrexGlass;
+	    energyLoss += 0.16 * dedx_Air;
+	}
+
+	return energyLoss;
+    }
+
+}