Author: [log in to unmask]
Date: Wed Jan 7 17:10:28 2015
New Revision: 1895
Log:
Add cluster energy and position calculation classes extracted from HPSEcalCluster and HPSEcalClusterIC.
Added:
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterEnergyCorrection.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPositionCorrection.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPropertyCalculator.java
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterEnergyCorrection.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterEnergyCorrection.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterEnergyCorrection.java Wed Jan 7 17:10:28 2015
@@ -0,0 +1,82 @@
+package org.hps.recon.ecal.cluster;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.base.BaseCluster;
+
+/**
+ * Cluster energy correction algorithm extracted from <code>HPSEcalClusterIC</code> class.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public final class ReconClusterEnergyCorrection {
+
+ // Variables for electron energy corrections.
+ static final double ELECTRON_ENERGY_A = -0.0027;
+ static final double ELECTRON_ENERGY_B = -0.06;
+ static final double ELECTRON_ENERGY_C = 0.95;
+
+ // Variables for positron energy corrections.
+ static final double POSITRON_ENERGY_A = -0.0096;
+ static final double POSITRON_ENERGY_B = -0.042;
+ static final double POSITRON_ENERGY_C = 0.94;
+
+ // Variables for photon energy corrections.
+ static final double PHOTON_ENERGY_A = 0.0015;
+ static final double PHOTON_ENERGY_B = -0.047;
+ static final double PHOTON_ENERGY_C = 0.94;
+
+ /**
+ * Calculate the corrected energy for the cluster.
+ * @param cluster The input cluster.
+ * @return The corrected energy.
+ */
+ public static double calculateCorrectedEnergy(Cluster cluster) {
+ double rawE = cluster.getEnergy();
+ return computeCorrectedEnergy(cluster.getParticleId(), rawE);
+ }
+
+ /**
+ * Calculate the corrected energy and set on the cluster.
+ * @param cluster The input cluster.
+ */
+ public static void setCorrectedEnergy(BaseCluster cluster) {
+ double correctedEnergy = calculateCorrectedEnergy(cluster);
+ cluster.setEnergy(correctedEnergy);
+ }
+
+ /**
+ * Calculates energy correction based on cluster raw energy and particle type as per
+ * <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param pdg Particle id as per PDG
+ * @param rawEnergy Raw Energy of the cluster (sum of hits with shared hit distribution)
+ * @return Corrected Energy
+ */
+ private static double computeCorrectedEnergy(int pdg, double rawEnergy) {
+ switch(pdg) {
+ case 11:
+ // electron
+ return computeCorrectedEnergy(rawEnergy, ELECTRON_ENERGY_A, ELECTRON_ENERGY_B, ELECTRON_ENERGY_C);
+ case -11:
+ // positron
+ return computeCorrectedEnergy(rawEnergy, POSITRON_ENERGY_A, POSITRON_ENERGY_B, POSITRON_ENERGY_C);
+ case 22:
+ // photon
+ return computeCorrectedEnergy(rawEnergy, PHOTON_ENERGY_A, PHOTON_ENERGY_B, PHOTON_ENERGY_C);
+ default:
+ // unknown
+ return rawEnergy;
+ }
+ }
+
+ /**
+ * Calculates the energy correction to a cluster given the variables from the fit as per
+ * <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param rawEnergy Raw energy of the cluster
+ * @param A,B,C from fitting in note
+ * @return Corrected Energy
+ */
+ private static double computeCorrectedEnergy(double rawEnergy, double varA, double varB, double varC){
+ double corrEnergy = rawEnergy / (varA * rawEnergy + varB / (Math.sqrt(rawEnergy)) + varC);
+ return corrEnergy;
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPositionCorrection.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPositionCorrection.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPositionCorrection.java Wed Jan 7 17:10:28 2015
@@ -0,0 +1,95 @@
+package org.hps.recon.ecal.cluster;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.base.BaseCluster;
+
+/**
+ * <p>
+ * Cluster position corrections taken from the <code>HPSEcalClusterIC</code> class.
+ * <p>
+ * This should be used before the energy is corrected on the Cluster.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public final class ReconClusterPositionCorrection {
+
+ // Variables for electron position corrections.
+ static final double ELECTRON_POS_A = 0.0066;
+ static final double ELECTRON_POS_B = -0.03;
+ static final double ELECTRON_POS_C = 0.028;
+ static final double ELECTRON_POS_D = -0.45;
+ static final double ELECTRON_POS_E = 0.465;
+
+ // Variables for positron position corrections.
+ static final double POSITRON_POS_A = 0.0072;
+ static final double POSITRON_POS_B = -0.031;
+ static final double POSITRON_POS_C = 0.007;
+ static final double POSITRON_POS_D = 0.342;
+ static final double POSITRON_POS_E = 0.108;
+
+ // Variables for photon position corrections.
+ static final double PHOTON_POS_A = 0.005;
+ static final double PHOTON_POS_B = -0.032;
+ static final double PHOTON_POS_C = 0.011;
+ static final double PHOTON_POS_D = -0.037;
+ static final double PHOTON_POS_E = 0.294;
+
+ public static double[] calculateCorrectedPosition(Cluster cluster) {
+ double clusterPosition[] = cluster.getPosition();
+ double correctedPosition = computeCorrectedPosition(cluster.getParticleId(), clusterPosition[0], cluster.getEnergy());
+ double[] position = new double[3];
+ position[0] = correctedPosition;
+ position[1] = clusterPosition[1];
+ position[2] = clusterPosition[2];
+ return position;
+ }
+
+ public static void setCorrectedPosition(BaseCluster cluster) {
+ cluster.setPosition(calculateCorrectedPosition(cluster));
+ }
+
+ /**
+ * Calculates position correction based on cluster raw energy, x calculated position,
+ * and particle type as per
+ * <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param pdg Particle id as per PDG
+ * @param xCl Calculated x centroid position of the cluster, uncorrected, at face
+ * @param rawEnergy Raw energy of the cluster (sum of hits with shared hit distribution)
+ * @return Corrected x position
+ */
+ private static double computeCorrectedPosition(int pdg, double xPos, double rawEnergy) {
+ double xCl = xPos / 10.0;//convert to mm
+ double xCorr;
+ switch(pdg) {
+ case 11: //Particle is electron
+ xCorr = positionCorrection(xCl, rawEnergy, ELECTRON_POS_A, ELECTRON_POS_B, ELECTRON_POS_C, ELECTRON_POS_D, ELECTRON_POS_E);
+ return xCorr * 10.0;
+ case -11:// Particle is positron
+ xCorr = positionCorrection(xCl, rawEnergy, POSITRON_POS_A, POSITRON_POS_B, POSITRON_POS_C, POSITRON_POS_D, POSITRON_POS_E);
+ return xCorr * 10.0;
+ case 22: // Particle is photon
+ xCorr = positionCorrection(xCl, rawEnergy, PHOTON_POS_A, PHOTON_POS_B, PHOTON_POS_C, PHOTON_POS_D, PHOTON_POS_E);
+ return xCorr * 10.0;
+ default: //Unknown
+ xCorr = xCl;
+ return xCorr * 10.0;
+ }
+ }
+
+ /**
+ * Calculates the position correction in cm using the raw energy and variables associated with the fit
+ * of the particle as described in
+ * <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param xCl
+ * @param rawEnergy
+ * @param varA
+ * @param varB
+ * @param varC
+ * @param varD
+ * @param varE
+ * @return
+ */
+ private static double positionCorrection(double xCl, double rawEnergy, double varA, double varB, double varC, double varD, double varE) {
+ return xCl - (varA / Math.sqrt(rawEnergy) + varB ) * xCl - (varC * rawEnergy + varD / Math.sqrt(rawEnergy) + varE);
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPropertyCalculator.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPropertyCalculator.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ReconClusterPropertyCalculator.java Wed Jan 7 17:10:28 2015
@@ -0,0 +1,319 @@
+package org.hps.recon.ecal.cluster;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.List;
+
+import org.hps.recon.ecal.ECalUtils;
+import org.lcsim.detector.IGeometryInfo;
+import org.lcsim.detector.solids.Trd;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.base.AbstractClusterPropertyCalculator;
+
+/**
+ * <p>
+ * Compute cluster properties for an HPS ECAL cluster.
+ * <p>
+ * This was extracted from the <code>HPSEcalCluster</code> class.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public final class ReconClusterPropertyCalculator extends AbstractClusterPropertyCalculator {
+
+ static final double eCriticalW = 800.0 * ECalUtils.MeV / (74 + 1);
+ static final double radLenW = 8.8; // mm
+
+ private double dmax;
+
+ /**
+ * Calculate the max distance parameter for use in the property calculation.
+ * @param energy The energy of the cluster.
+ */
+ private void calculateMaxDistance(double energy, boolean isElectron) {
+ double y = energy / eCriticalW;
+ double Cj = isElectron ? -0.5 : 0.5;
+ double tmax = Math.log(y) + Cj; //Maximum of dE/dt profile in units of rad. len.
+ dmax = tmax * radLenW; //mm
+ }
+
+ /**
+ * <p>
+ * Perform the cluster property calculations, which sets this object's state.
+ * <p>
+ * Copied and modified from
+ * {@link org.lcsim.event.base.TensorClusterPropertyCalculator#calculateProperties(List)}.
+ */
+ public void calculateProperties(Cluster cluster) {
+
+ calculateMaxDistance(cluster.getEnergy(), cluster.getParticleId() == 11 || cluster.getParticleId() == 0);
+
+ List<CalorimeterHit> hits = cluster.getCalorimeterHits();
+
+ double[] mm_NE = new double[3];
+ double[] mm_CE = new double[3];
+ double[][] 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);
+
+ //Find position at shower max
+ IGeometryInfo geom = hit.getDetectorElement().getGeometry();
+ double[] pos = geom.transformLocalToGlobal(VecOp.add(geom.transformGlobalToLocal(geom.getPosition()), (Hep3Vector) new BasicHep3Vector(0, 0, dmax - 1 * ((Trd) geom.getLogicalVolume().getSolid()).getZHalfLength()))).v();
+
+ //pos = global_pos_tmax.v();
+ double E = hit.getCorrectedEnergy();
+ 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 eigenvalue EV, the axis is (nx, ny, nz) where:
+ // (Exx - EV)nx + (Exy)ny + (Exz)nz = 0
+ // (Eyx)nx + (Eyy - EV)ny + (Eyz)nz = 0
+ // (Ezx)nx + (Ezy)ny + (Ezz - EV)nz = 0
+ // Setting nx = 1, we have:
+ // (Exx - EV) + (Exy)ny + (Exz)nz = 0
+ // (Eyx) + (Eyy - EV)ny + (Eyz)nz = 0
+ // (Ezx) + (Ezy)ny + (Ezz - EV)nz = 0
+ // and so
+ // (Exy)ny = EV - Exx - (Exz)nz => ny = (EV - Exx - Exz*nz)/Exy
+ // What if Exy = 0? Then provided Eyz is non-zero we can write:
+ // (Ezx) + (Ezy)ny + (Ezz - EV)nz = 0
+ // ny = (Exz - (Ezz-EV)*nz)/Eyz
+ // What if Exy = 0 and Eyz = 0 but (Eyy - EV) is non-zero?
+ // (Eyy - EV)ny + (Eyz)nz = 0
+ // ny = -(Eyz*nz)/(Eyy-EV)
+
+ // In the pathological case where Exz = Eyz = Ezz = 0:
+ // (Exx - EV)nx + (Exy)ny = 0 => ny/nx = -(Exx-EV)/Exy
+ // (Eyx)nx + (Eyy - EV)ny = 0 => ny/nx = -Eyx/(Eyy-EV)
+ // (EV)nz = 0
+ // so
+ // -ny/nx = (EV-Exx)/Exy = Eyx/(EV-Eyy)
+ // But watch out for order! Recalculate eigenvalues for this pathological case.
+ // (EV-Exx)(EV-Eyy) = Eyx*Exy
+ // EV^2 - EV(Exx+Eyy) + Exx*Eyy - Eyx*Exy = 0
+ //
+ // In another pathological case, Exz = Exy = 0:
+ // (Exx - EV)nx = 0
+ // (Eyy - EV)ny + (Eyz)nz = 0 => ny/nz = -(Eyz)/(Eyy-EV)
+ // (Ezy)ny + (Ezz - EV)nz = 0 => ny/nz = -(Ezz-EV)/(Ezy)
+ // so we cannot set nx = 1. Instead, write:
+ // -ny/nz = (Eyz)/(Eyy-EV) = (Ezz-EV)/(Ezy)
+ // Then
+ // (Eyz)(Ezy) = (Eyy-EV)(Ezz-EV)
+ // (Eyz)^2 = (Eyy)(Ezz) - (Eyy)(EV) - (Ezz)(EV) + (EV)^2
+ // EV^2 - EV(Eyy+Ezz) + Eyy*Ezz - Eyz*Eyz = 0
+ // Handle pathological case
+ if (Exz == 0.0 && Eyz == 0.0) {
+ // Recompute eigenvectors.
+ EV[0] = 0.5 * (Exx + Eyy) + 0.5 * Math.sqrt((Exx + Eyy) * (Exx + Eyy) + 4.0 * Exy * Exy);
+ EV[1] = 0.5 * (Exx + Eyy) - 0.5 * Math.sqrt((Exx + Eyy) * (Exx + Eyy) + 4.0 * Exy * Exy);
+ EV[2] = 0.0;
+ for (int i = 0; i < 2; i++) {
+ double nx_over_ny = Exy / (Exx - EV[i]);
+ double nx_unnormalized = nx_over_ny;
+ double ny_unnormalized = 1.0;
+ double norm = Math.sqrt(nx_unnormalized * nx_unnormalized + ny_unnormalized * ny_unnormalized);
+ mm_PA[i][0] = ny_unnormalized / norm;
+ mm_PA[i][1] = nx_unnormalized / norm;
+ mm_PA[i][2] = 0.0;
+ }
+ // ... and now set third eigenvector to the z direction:
+ mm_PA[2][0] = 0.0;
+ mm_PA[2][1] = 0.0;
+ mm_PA[2][2] = 1.0;
+ } else if (Exz == 0.0 && Exy == 0.0) {
+ // Another pathological case
+ EV[0] = 0.5 * (Eyy + Ezz) + 0.5 * Math.sqrt((Eyy + Ezz) * (Eyy + Ezz) + 4.0 * Eyz * Eyz);
+ EV[1] = 0.5 * (Eyy + Ezz) - 0.5 * Math.sqrt((Eyy + Ezz) * (Eyy + Ezz) + 4.0 * Eyz * Eyz);
+ EV[2] = 0.0;
+ for (int i = 0; i < 2; i++) {
+ double ny_over_nz = Eyz / (Eyy - EV[i]);
+ double ny_unnormalized = ny_over_nz;
+ double nz_unnormalized = 1.0;
+ double norm = Math.sqrt(ny_unnormalized * ny_unnormalized + nz_unnormalized * nz_unnormalized);
+ mm_PA[i][0] = nz_unnormalized / norm;
+ mm_PA[i][1] = ny_unnormalized / norm;
+ mm_PA[i][2] = 0.0;
+ }
+ mm_PA[2][0] = 0.0;
+ mm_PA[2][1] = 0.0;
+ mm_PA[2][2] = 1.0;
+ } else {
+ 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;
+ if (Exy == 0.0) {
+ // Recompute
+ if (Eyz != 0.0) {
+ // ny = (Exz - (Ezz-EV)*nz)/Eyz
+ C[1] = (Exz - (Ezz - EV[i]) * C[2]) / Eyz;
+ } else {
+ // ny = -(Eyz*nz)/(Eyy-EV)
+ C[1] = -(Eyz * C[2]) / (Eyy - EV[i]);
+ }
+ }
+ 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];
+ }
+
+ // Compute theta angle from origin to cluster center.
+ itheta = Math.atan(position[1] / position[2]);
+
+ // Compute phi angle from origin to cluster center.
+ iphi = Math.atan(position[0] / position[2]);
+ }
+}
