diff -u -r1.32 -r1.33
--- ReclusterDTreeDriver.java	17 Jul 2008 01:03:32 -0000	1.32
+++ ReclusterDTreeDriver.java	20 Jul 2008 22:50:03 -0000	1.33
@@ -34,7 +34,7 @@

   * in this package, which uses the implementation in
   * org.lcsim.recon.cluster.directedtree developed by NIU).
   *

-  * @version $Id: ReclusterDTreeDriver.java,v 1.32 2008/07/17 01:03:32 mcharles Exp $

+  * @version $Id: ReclusterDTreeDriver.java,v 1.33 2008/07/20 22:50:03 mcharles Exp $

   * @author Mat Charles <[log in to unmask]>
   */
 

@@ -127,6 +127,12 @@

 	super.debugProcess(event);
 	m_event = event;
 

+	// Steve's pre-shower MIP-finder
+	{
+	    SteveMipWrapper tmpWrapper = new SteveMipWrapper();
+	    tmpWrapper.process(event);
+	}
+

 	if (m_oldMipFinderCrossesTrees) {
 	    m_allowComponentsToStraddleLargeClusters = true;
 	}

@@ -579,8 +585,6 @@

 	    sharedLeftoverHitClustersHCAL.rebuildHints();
 	    allSharedClusters.add(sharedLeftoverHitClustersHCAL);
 	} else if (steveMipShareForHaloHCAL) {

-	    SteveMipWrapper tmpWrapper = new SteveMipWrapper();
-	    tmpWrapper.process(event);

 	    SteveMIPReassignmentAlgorithm tmpMipAlg = new SteveMIPReassignmentAlgorithm(event, 1.0);
 	    DownstreamTrackClusterSharingAlgorithm coneSharingAlgHCAL  = new DownstreamTrackClusterSharingAlgorithm(clustersMatchedToTweakedTracks, tmpMipAlg); // TEST
 	    SharedClusterGroup sharedLeftoverHitClustersHCAL = new SharedClusterGroup(leftoverHitClustersToShareHCAL, coneSharingAlgHCAL); // TEST

diff -u -r1.3 -r1.4
--- ChargedHadronClusterEnergyCalculator.java	27 Mar 2008 23:34:07 -0000	1.3
+++ ChargedHadronClusterEnergyCalculator.java	20 Jul 2008 22:50:03 -0000	1.4
@@ -1,6 +1,8 @@

 package org.lcsim.recon.pfa.structural;
 
 import java.util.*;

+import hep.physics.vec.*;
+

 import org.lcsim.recon.cluster.util.*;
 import org.lcsim.event.Cluster;
 import org.lcsim.event.CalorimeterHit;

@@ -28,11 +30,9 @@

   * Since the number of non-MIP hits could be small or zero, the neutral calibration should not have
   * a large zero offset.
   *

-  * Still to fix:
-  *   1) Increase energy loss for tracks that pass through the layer at an angle (currently
-  *      it assumes all tracks are perpendicular to the material planes).

+  * @author Mat Charles <[log in to unmask]>

   *

-  * @version $Id: ChargedHadronClusterEnergyCalculator.java,v 1.3 2008/03/27 23:34:07 mcharles Exp $

+  * @version $Id: ChargedHadronClusterEnergyCalculator.java,v 1.4 2008/07/20 22:50:03 mcharles Exp $

   */
 
 public class ChargedHadronClusterEnergyCalculator extends Driver implements ClusterEnergyCalculator

@@ -40,46 +40,71 @@

     ClusterEnergyCalculator m_neutralCalib = null;
     String m_mipListName = null;
     EventHeader m_event = null;

+    Map<Cluster, Hep3Vector> m_mipDirectionCache;
+

     public ChargedHadronClusterEnergyCalculator(String mipListName, ClusterEnergyCalculator neutralCalib) {
 	m_neutralCalib = neutralCalib;
 	m_mipListName = mipListName;
     }
     public void process(EventHeader event) {
 	m_event = event;

+	m_mipDirectionCache = new HashMap<Cluster, Hep3Vector>(); // Flush

     }
     public double getEnergy(Cluster c){
 	// Grab MIP list from the event:

-	List<Cluster> mips = m_event.get(Cluster.class, m_mipListName);

+	List<Cluster> mips = null;
+	if (m_event.hasCollection(Cluster.class, m_mipListName)) {
+	    mips = m_event.get(Cluster.class, m_mipListName);
+	} else {
+	    System.out.println("ERROR: Event has no MIP list (name='"+m_mipListName+"'). Unless this is a single-particle event, fix this error!");
+	    mips = new Vector<Cluster>(); // empty
+	}

 	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();

+	BasicCluster neutralHadronCluster = 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);

+		neutralHadronCluster.addHit(hit);

 	    } else {
 		// From a MIP segment
 		mipHitsInCluster.add(hit);

-		mipEnergySum += getMIPEnergy(hit);

 	    }
 	}

-	double neutralHadronEnergy = m_neutralCalib.getEnergy(tmpCopy);
-	return mipEnergySum+neutralHadronEnergy;
-    }
-    private double getMIPEnergy(Collection<CalorimeterHit> mipHits) {
-	double eSum = 0.0;
-	for (CalorimeterHit hit : mipHits) {
-	    eSum += getMIPEnergy(hit);

+
+
+	BasicCluster mipHitsClus = new BasicCluster();
+	for (CalorimeterHit hit : mipHitsInCluster) { mipHitsClus.addHit(hit); }
+	Hep3Vector direction = getMIPDirectionFromOrigin(mipHitsClus);
+
+	boolean useSimpleDirection = true;
+	double mipEnergySum = 0.0;
+	for (CalorimeterHit hit : mipHitsInCluster) {
+	    double hitWeight = 1.0;
+	    if (hit instanceof FuzzyCalorimeterHit) {
+		hitWeight = ((FuzzyCalorimeterHit)(hit)).getWeight();
+	    }
+	    double hitEnergyAtNormalIncidence = getMIPEnergyAtNormalIncidence(hit);
+	    Hep3Vector estimatedDirection = direction;
+	    if (!useSimpleDirection) {
+		estimatedDirection = getMIPDirectionFromSubClusters(hit, mips);
+	    }
+	    Hep3Vector estimatedNormalIncidence = getNormal(hit);
+	    double cosTheta = VecOp.dot(VecOp.unit(estimatedDirection), VecOp.unit(estimatedNormalIncidence));
+	    double angularCorrection = getAngularCorrection(cosTheta);
+	    double hitEnergy = hitEnergyAtNormalIncidence * angularCorrection;
+	    mipEnergySum += hitWeight * hitEnergy;

 	}

-	return eSum;

+
+	double neutralHadronEnergy = m_neutralCalib.getEnergy(neutralHadronCluster);
+	return mipEnergySum+neutralHadronEnergy;

     }
 

-    private double getMIPEnergy(CalorimeterHit hit) {

+    private double getMIPEnergyAtNormalIncidence(CalorimeterHit hit) {

 	// First, grab hit info (ID, layer)
 	org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
 	id.setID(hit.getCellID());

@@ -139,8 +164,156 @@

 		throw new AssertionError("Don't know how to get layer thickness for geometry type "+subLayerSolid.getClass().getName());
 	    }
 	}

-

+	

 	return estimatedEnergyLoss;
     }
 

+    Hep3Vector getNormal(CalorimeterHit hit) {
+	org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+	if (subdet.isEndcap()) {
+	    return new BasicHep3Vector(0.0, 0.0, 1.0);
+	} else {
+	    // For barrel cylinder, normal incidence is given by the
+	    // vector from the z-axis to the cell position. The vector
+	    // has no z-component. The x and y components are determined
+	    // by the x and y co-ordinates of the cell position.
+	    Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+	    double x = hitPosition.x();
+	    double y = hitPosition.y();
+	    Hep3Vector tmpVector = new BasicHep3Vector(x, y, 0.0);
+	    return VecOp.unit(tmpVector);
+	}
+    }
+
+    // Now, all of the above assumed normal incidence. In practice, the particle
+    // was probably travelling at some angle theta. Then the actual distance
+    // covered is increased by a factor of 1/cos(theta) where cos(theta) is the
+    // angle between normal incidence and the track direction
+
+    double getAngularCorrection(double cosTheta) {
+	// Make sure we're working with a positive number:
+	cosTheta = Math.abs(cosTheta);
+
+	if (cosTheta == 0.0 || Double.isNaN(cosTheta)) {
+	    // Something broke. This shouldn't happen.
+	    // Panic, then recover semi-gracefully by assuming normal incidence
+	    return 1.0;
+	}
+	boolean useRonScaling = false;
+	double scaleFactor = getSimpleScaleFactor(cosTheta);
+	if (useRonScaling) {
+	    scaleFactor = getRonScaleFactor(cosTheta);
+	}
+	if (scaleFactor < 1.0 || Double.isNaN(scaleFactor) || scaleFactor > 100.0) {
+	    // Something broke -- treat it as normal incidence
+	    return 1.0;
+	} else {
+	    // All OK
+	    return scaleFactor;
+	}
+    }
+
+    double getSimpleScaleFactor(double cosTheta) {
+	return 1.0 / cosTheta; // very simple scaling
+    }
+
+    double getRonScaleFactor(double cosTheta) {
+	// Scale according to Ron's formula:
+	//    1 / (1. + alpha*(1./st - 1.));
+	// where st = sin(theta) or cos(theta)
+	// and alpha is a constant (hard-code to -0.23 for now... can read from hadronCalibration/[...] in properties?
+	double alpha = -0.23;
+	double scaleFactor = 1.0 / (1.0 + alpha*(1.0/cosTheta - 1.0));
+	return scaleFactor;
+    }
+
+    // Now estimate the direction the MIP was travelling in so we can determine the angle of incidence.
+    //
+    // My first try at this involved looking at the MIP hits and estimating
+    // a direction based on the principal axes of a cluster made from those hits.
+    // However, this can go badly wrong:
+    //   a) if more than one MIP contributes to the shower -- then the MIP hits
+    //      will form a weirdly-shaped blob and the principal axes mean nothing;
+    //   b) if the principal axis calculation goes wrong (can happen in a few
+    //      pathological cases).
+    // It can also go wrong if there are few/no MIP hits, but that is easy to
+    // spot.
+    //
+    // Because problem (a) is serious and hard to detect, we'll go with a very
+    // simple approximation: the MIP is assumed to be coming from the IP in a
+    // straight line. This is not particularly good but it is at least robust.
+
+
+    Hep3Vector getMIPDirectionFromOrigin(Cluster clus) {
+	if (clus.getCalorimeterHits().size()==0) {
+	    // Direction not defined
+	    return null;
+	} else {
+	    // Assume it flew in a straight line from the origin.
+	    Hep3Vector pos = new BasicHep3Vector(clus.getPosition());
+	    return VecOp.unit(pos);
+	}
+    }
+
+    // As an alternative, we can treat each MIP separately (for the case of multiple
+    // MIP clusters). Then the direction of each MIP will be better-defined and we
+    // can correct for hits in the different MIPs according to the direction of
+    // that particular MIP.
+    //
+    // Watch out! Depending on the MIP-finding algorithm, the same hit might appear
+    // in more than one MIP! We must take care not to double-count.
+
+    Hep3Vector getMIPDirectionFromSubClusters(CalorimeterHit hit, List<Cluster> mips) {
+	List<Cluster> matchedMips = new Vector<Cluster>();
+	for (Cluster mip : mips) {
+	    if (mip.getCalorimeterHits().contains(hit)) {
+		matchedMips.add(mip);
+	    }
+	}
+	// Make sure we look at one and only one MIP:
+	Cluster uniqueMatchedMip = null;
+	if (matchedMips.size() < 1) { 
+	    throw new AssertionError("Book-keeping error! Failed to match hit to any mips.");
+	} else if (matchedMips.size() == 1) {
+	    // Unambiguous match
+	    uniqueMatchedMip = matchedMips.get(0);
+	} else {
+	    // Matched to >1 MIP. Use the one with the most hits (since this will have the
+	    // best direction estimate):
+	    for (Cluster mip : matchedMips) {
+		if (uniqueMatchedMip==null || uniqueMatchedMip.getCalorimeterHits().size()<mip.getCalorimeterHits().size()) {
+		    // Found a bigger one
+		    uniqueMatchedMip = mip;
+		}
+	    }
+	}
+
+	// Now find the direction of the MIP:
+	return getMIPDirectionFromPointing(uniqueMatchedMip);
+    }
+
+    Hep3Vector getMIPDirectionFromPointing(Cluster clus) {
+	if (clus.getCalorimeterHits().size()<4) {
+	    // Not enough hits to define a direction
+	    // Fall back to simple method
+	    return getMIPDirectionFromOrigin(clus);
+	} else {
+	    // Enough hits to define a direction by fitting as a line.
+	    Hep3Vector cachedDirection = m_mipDirectionCache.get(clus);
+	    if (cachedDirection != null) {
+		return cachedDirection;
+	    } else {
+		// Make a clone cluster that we can safely calculate properties on:
+		TensorClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
+		BasicCluster copyOfClus = new BasicCluster();
+		for (CalorimeterHit hit : clus.getCalorimeterHits()) { copyOfClus.addHit(hit); }
+		copyOfClus.setPropertyCalculator(calc);
+		copyOfClus.calculateProperties();
+		double[][] axes = calc.getPrincipleAxis();
+		Hep3Vector direction = VecOp.unit(new BasicHep3Vector(axes[0][0], axes[0][1], axes[0][2]));
+		m_mipDirectionCache.put(clus, direction);
+		return direction;
+	    }
+	}
+    }

 }