diff -N TensorClusterPropertyCalculator.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ TensorClusterPropertyCalculator.java	30 Jun 2005 20:55:23 -0000	1.1
@@ -0,0 +1,280 @@

+package org.lcsim.recon.cluster.util;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+
+/**
+ * Implementation of ClusterPropertyCalculator Interface using
+ * inertia Tensor calculation to derive parameters 
+ * 
+ * @author cassell
+ */
+public class TensorClusterPropertyCalculator implements ClusterPropertyCalculator
+{
+	protected double[] position;
+	protected double[] positionError;
+	protected double iphi;
+	protected double itheta;
+	protected double[] directionError;
+	protected double[] shapeParameters;
+	protected double[] mm_NE;
+                             protected double[] mm_CE;
+                             protected double[][] mm_PA;
+
+	/** 
+	 * Calculate the cluster properties from the hits
+	 */
+	public TensorClusterPropertyCalculator()
+	{
+		position = new double[3];
+		positionError = new double[6];
+		iphi = 0.;
+		itheta = 0.;
+		directionError = new double[6];
+		shapeParameters = new double[6];
+	}
+	public void calculateProperties(List<CalorimeterHit> hits)
+	{
+		mm_NE = new double[3];
+		mm_CE = new double[3];
+		mm_PA = new double[3][3]; 
+		for(int i=0;i<3;++i)
+		{
+			mm_NE[i] = 0.;
+			mm_CE[i] = 0.;
+			for(int j=0;j<3;++j){mm_PA[i][j]= 0.;}
+		}
+		double Etot = 0.0;
+		double Exx = 0.0;
+		double Eyy = 0.0;
+		double Ezz = 0.0;
+		double Exy = 0.0;
+		double Eyz = 0.0;
+		double Exz = 0.0;
+		double CEx = 0.0;
+		double CEy = 0.0;
+		double CEz = 0.0;
+		double CEr = 0.0;
+		double E1 = 0.0;
+		double E2 = 0.0;
+		double E3 = 0.0;
+		double NE1 = 0.0;
+		double NE2 = 0.0;
+		double NE3 = 0.0;
+		double Tr = 0.0;
+		double M = 0.0;
+		double Det = 0.0;
+		int nhits = hits.size();
+		for(int i=0;i<hits.size();i++)
+		{
+			CalorimeterHit hit = hits.get(i);
+			double[] pos = hit.getPosition();
+		//	Subdetector d = hit.getSubdetector();
+		//	double sampling_fraction = d.getSamplingFraction();
+			double sampling_fraction = 1.;
+			double E = hit.getEnergy()/sampling_fraction;
+			Etot += E;
+			CEx += E*pos[0];
+			CEy += E*pos[1];
+			CEz += E*pos[2];
+			Exx += E*(pos[1]*pos[1] + pos[2]*pos[2]);
+			Eyy += E*(pos[0]*pos[0] + pos[2]*pos[2]);
+			Ezz += E*(pos[1]*pos[1] + pos[0]*pos[0]);
+			Exy -= E*pos[0]*pos[1];
+			Eyz -= E*pos[1]*pos[2];
+			Exz -= E*pos[0]*pos[2];
+		}
+		CEx = CEx/Etot;
+		CEy = CEy/Etot;
+		CEz = CEz/Etot;
+		double CErSq = CEx*CEx + CEy*CEy + CEz*CEz;
+		CEr = Math.sqrt(CErSq);
+		// now go to center of energy coords.
+		if (nhits > 3 )
+		{
+			Exx = Exx - Etot*(CErSq - CEx*CEx);
+			Eyy = Eyy - Etot*(CErSq - CEy*CEy);
+			Ezz = Ezz - Etot*(CErSq - CEz*CEz);
+			Exy = Exy + Etot*CEx*CEy;
+			Eyz = Eyz + Etot*CEy*CEz;
+			Exz = Exz + Etot*CEz*CEx;
+
+	//
+			Tr = Exx + Eyy + Ezz;
+			double Dxx = Eyy*Ezz - Eyz*Eyz;
+			double Dyy = Ezz*Exx - Exz*Exz;
+			double Dzz = Exx*Eyy - Exy*Exy;
+			M = Dxx + Dyy + Dzz;
+			double Dxy = Exy*Ezz - Exz*Eyz;
+			double Dxz = Exy*Eyz - Exz*Eyy;
+			Det = Exx*Dxx - Exy*Dxy + Exz*Dxz;
+			double xt = Tr*Tr - 3*M;
+			double sqrtxt = Math.sqrt(xt);
+	// eqn to solve for eigenvalues is x**3 - x**2*Tr + x*M - Det = 0
+	// crosses y axis at -Det and inflection points are 
+				
+			double mE1 = 0.;
+			double mE2 = 0.;
+			double mE3 = 0.;
+			double a = (3*M - Tr*Tr)/3.;
+			double b = (-2*Tr*Tr*Tr + 9*Tr*M -27*Det)/27.;
+			double test = b*b/4. + a*a*a/27.;
+			if(test >= 0.01)
+			{
+				System.out.println("AbstractCluster: Only 1 real root!!!");
+				System.out.println("  nhits = " + nhits + "\n");
+				System.out.println(" a,b,test = " + a + " " + b + " " + test + "\n");
+			}
+			else
+			{
+				double temp;
+				if(test >= 0.)temp = 1.;
+				else temp = Math.sqrt(b*b*27./(-a*a*a*4.));
+				if(b > 0.)temp = -temp;
+				double phi = Math.acos(temp);
+				double temp1 = 2.*Math.sqrt(-a/3.);
+				mE1 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3.);
+				mE2 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3. + 2.*Math.PI/3.);
+				mE3 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3. + 4.*Math.PI/3.);
+			}
+			if(mE1 < mE2)
+			{
+				if(mE1 < mE3)
+				{
+					E1 = mE1;
+					if(mE2 < mE3)
+					{
+						E2 = mE2;
+						E3 = mE3;
+					}
+					else
+					{
+						E2 = mE3;
+						E3 = mE2;
+					}
+				}
+				else
+				{
+					E1 = mE3;
+					E2 = mE1;
+					E3 = mE2;
+				}
+			}
+			else
+			{
+				if(mE2 < mE3)
+				{
+					E1 = mE2;
+					if(mE1 < mE3)
+					{
+						E2 = mE1;
+						E3 = mE3;
+					}
+					else
+					{
+						E2 = mE3;
+						E3 = mE1;
+					}
+				}
+				else
+				{
+					E1 = mE3;
+					E2 = mE2;
+					E3 = mE1;
+				}
+			}
+			
+			NE1 = E1/Etot;
+			NE2 = E2/Etot;
+			NE3 = E3/Etot;
+			double[] EV = new double[3];
+			EV[0] = E1;
+			EV[1] = E2;
+			EV[2] = E3;
+	// Now calculate principal axes
+			for( int i = 0 ; i < 3 ; i++ )
+			{
+				double[] C = new double[3];
+				C[0] = 1.0;
+				C[2] = (Exy*Exy + (Eyy - EV[i])*(EV[i] - Exx))/
+				((Eyy - EV[i])*Exz - Eyz*Exy);
+				C[1] = (EV[i] - Exx - Exz*C[2])/Exy;
+				double norm = Math.sqrt(C[0]*C[0] + C[1]*C[1] + C[2]*C[2]);
+				mm_PA[i][0] = C[0]/norm;
+				mm_PA[i][1] = C[1]/norm;
+				mm_PA[i][2] = C[2]/norm;
+			}
+		}
+		mm_NE[0] = NE1;
+		mm_NE[1] = NE2;
+		mm_NE[2] = NE3;
+		mm_CE[0] = CEx;
+		mm_CE[1] = CEy;
+		mm_CE[2] = CEz;
+		for(int i=0;i<6;i++)
+		{
+			positionError[i] = 0.;
+			directionError[i] = 0.;
+ 			shapeParameters[i] = 0.;
+		}
+		for(int i=0;i<3;i++)
+		{
+			position[i] = mm_CE[i];
+			shapeParameters[i] = mm_NE[i];
+		}
+		if(nhits > 3)
+		{
+			double dr = Math.sqrt(  (position[0]+mm_PA[0][0])*(position[0]+mm_PA[0][0]) +
+									(position[1]+mm_PA[0][1])*(position[1]+mm_PA[0][1]) +
+									(position[2]+mm_PA[0][2])*(position[2]+mm_PA[0][2]) ) -
+					    Math.sqrt(	(position[0])*(position[0]) +
+									(position[1])*(position[1]) +
+									(position[2])*(position[2]) ) ;
+			double sign = 1.;
+			if(dr < 0.)sign = -1.;
+			itheta = Math.acos(sign*mm_PA[0][2]);
+			iphi = Math.atan2(sign*mm_PA[0][1],sign*mm_PA[0][0]);
+		}
+		else
+		{
+			itheta = 999.;
+			iphi = 999.;
+		}
+	}
+	public double[] getPosition()
+	{
+		return position;
+	}
+	public double[] getPositionError()
+	{
+		return positionError;
+	}
+	public double getIPhi()
+	{
+		return iphi;
+	}
+	public double getITheta()
+	{
+		return itheta;
+	}
+	public double[] getDirectionError()
+	{
+		return directionError;
+	}
+	public double[] getShapeParameters()
+	{
+		return shapeParameters;
+	}
+	public double[] getNormalizedEigenvalues()
+	{
+		return mm_NE;
+	}
+	public double[][] getPrincipleAxis()
+	{
+		return mm_PA;
+	}
+		
+}