Author: [log in to unmask]
Date: Wed Jan 7 17:03:40 2015
New Revision: 3479
Log:
Overhaul of BaseCluster class and related changes.
Modified:
projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/Cluster.java
projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/AbstractClusterPropertyCalculator.java
projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/BaseCluster.java
projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/ClusterPropertyCalculator.java
projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/TensorClusterPropertyCalculator.java
Modified: projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/Cluster.java
=============================================================================
--- projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/Cluster.java (original)
+++ projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/Cluster.java Wed Jan 7 17:03:40 2015
@@ -1,10 +1,8 @@
package org.lcsim.event;
import java.util.List;
-/** The LCIO cluster.
- *
- * @author gaede
- * @version $Id: Cluster.java,v 1.6 2011/11/09 22:50:28 jeremy Exp $
+/**
+ * The LCIO cluster interface.
*/
public interface Cluster
@@ -14,6 +12,8 @@
* elsewhere, e.g. in the run header. Bits 16-31 are used internally.
*/
public int getType();
+
+ public int getParticleId();
/** Energy of the cluster.
*/
Modified: projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/AbstractClusterPropertyCalculator.java
=============================================================================
--- projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/AbstractClusterPropertyCalculator.java (original)
+++ projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/AbstractClusterPropertyCalculator.java Wed Jan 7 17:03:40 2015
@@ -3,6 +3,7 @@
import java.util.List;
import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
/**
* An abstract implementation of {@link ClusterPropertyCalculator}.
@@ -11,15 +12,30 @@
*/
public abstract class AbstractClusterPropertyCalculator implements ClusterPropertyCalculator {
- protected double[] position;
- protected double[] positionError;
+ protected double[] position = new double[3];
+ protected double[] positionError = new double[6];
protected double iphi;
protected double itheta;
- protected double[] directionError;
- protected double[] shapeParameters;
+ protected double[] directionError = new double[6];
+ protected double[] shapeParameters = new double[6];
+
+ @Override
+ public void calculateProperties(List<CalorimeterHit> hits) {
+ }
@Override
- public abstract void calculateProperties(List<CalorimeterHit> hits);
+ public void calculateProperties(Cluster cluster) {
+ calculateProperties(cluster.getCalorimeterHits());
+ }
+
+ protected void reset() {
+ position = new double[3];
+ positionError = new double[6];
+ iphi = 0;
+ itheta = 0;
+ directionError = new double[6];
+ shapeParameters = new double[6];
+ }
@Override
public double[] getPosition() {
@@ -42,7 +58,7 @@
}
@Override
- public double[] getDirectionError() {
+ public double[] getDirectionError() {
return directionError;
}
Modified: projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/BaseCluster.java
=============================================================================
--- projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/BaseCluster.java (original)
+++ projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/BaseCluster.java Wed Jan 7 17:03:40 2015
@@ -7,385 +7,525 @@
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.Cluster;
-import org.lcsim.event.MCParticle;
-import org.lcsim.event.SimCalorimeterHit;
-import org.lcsim.geometry.IDDecoder;
-import org.lcsim.geometry.compact.Subdetector;
-
/**
- * Default implementation of Cluster Interface for Simulation
- *
- * @author cassell
+ * <p>
+ * This is a concrete implementation of the {@link org.lcsim.event.Cluster} LCIO interface.
+ * <p>
+ * This version is an overhaul of the previous base class, with the following changes:
+ * <ul>
+ * <li>Added several constructors with argument lists, including one that is fully qualified and another
+ * that takes a list of hits only.</li>
+ * <li>Removed the prior implementation of subdetector energies in favor of a simple set method.
+ * (This part of the API is basically unused anyways.)</li>
+ * <li>Added set methods for all class variables, where they were missing.</li>
+ * <li>Added access to particle ID based on the official API.</li>
+ * <li>Added a few utility methods for adding lists of hits and clusters.</li>
+ * <li>Added a method for setting a hit contribution that is different from the corrected energy.</li>
+ * <li>Simplified the {@link #calculateProperties()} method so it doesn't do a bunch of array copies.</li>
+ * </ul>
+ *
+ * @see org.lcsim.event.Cluster
+ * @see org.lcsim.event.CalorimeterHit
+ * @see org.lcsim.event.base.ClusterPropertyCalculator
+ * @see org.lcsim.event.base.TensorClusterPropertyCalculator
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @author Norman Graf <[log in to unmask]>
+ *
*/
-public class BaseCluster implements Cluster
-{
+public class BaseCluster implements Cluster {
+
protected List<CalorimeterHit> hits = new ArrayList<CalorimeterHit>();
+ protected List<Double> hitContributions = new ArrayList<Double>();
protected List<Cluster> clusters = new ArrayList<Cluster>();
- protected List<Subdetector> detectors = new ArrayList<Subdetector>();
-// protected double raw_energy;
- protected double corrected_energy;
+
+ protected double energy;
protected double energyError;
- protected double[] subdetector_raw_energies = new double[10];
- protected double[] subdetector_corrected_energies = new double[10];
- protected List<Double> hit_energies = new ArrayList<Double>();
- protected ClusterPropertyCalculator cluster_property_calculator = new TensorClusterPropertyCalculator();
- protected boolean needsPropertyCalculation = true;
- protected double[] position;
- protected double[] positionError;
+
+ protected double[] position = new double[3];
+ protected double[] positionError = new double[6];
+
protected double iphi;
protected double itheta;
- protected double[] directionError;
+
+ protected double[] directionError = new double[6];
protected double[] shapeParameters;
-
- public BaseCluster()
- {
-// raw_energy = 0.;
- corrected_energy = 0.;
- for(int i=0;i<10;i++)
- {
- subdetector_raw_energies[i] = 0.;
- subdetector_corrected_energies[i] = 0.;
- }
- position = new double[3];
- positionError = new double[6];
- directionError = new double[6];
- }
-
- /**
- * Add a CalorimeterHit to the cluster
- */
- public void addHit(CalorimeterHit hit)
- {
- hits.add(hit);
- double hre = hit.getRawEnergy();
-// raw_energy += hre;
- double hce = hit.getCorrectedEnergy();
- IDDecoder idc = hit.getIDDecoder();
- idc.setID(hit.getCellID());
- // Replace call to getValue("system").
- // -- 2007-04-12 Jeremy McCormick
- int detector_index = idc.getSystemID();
- if((detector_index > 9 )|(detector_index < 0))detector_index = 0;
- corrected_energy += hce;
- hit_energies.add(hce);
- subdetector_raw_energies[detector_index] += hre;
- subdetector_corrected_energies[detector_index] += hce;
+
+ protected double[] subdetectorEnergies = new double[0];
+
+ protected int type;
+ protected int pid;
+
+ protected ClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
+ protected boolean needsPropertyCalculation = true;
+
+ /**
+ * The no argument constructor.
+ */
+ public BaseCluster() {
+ }
+
+ /**
+ * Almost fully qualified constructor, if the cluster's properties are already calculated.
+ * The energy given here will override the value calculated from the hits. If this is not
+ * desired, then another constructor should be used instead. This constructor does not
+ * allow setting hit contributions that are different from the hit energies.
+ * @param hits The list of hits.
+ * @param energy The total energy.
+ * @param energyError The energy error.
+ * @param position The position.
+ * @param positionError The position error.
+ * @param iphi The intrinsic phi.
+ * @param itheta The intrinsic theta.
+ * @param directionError The direction error.
+ * @param shapeParameters The shape parameters.
+ */
+ public BaseCluster(
+ List<CalorimeterHit> hits,
+ double energy,
+ double energyError,
+ double[] position,
+ double[] positionError,
+ double iphi,
+ double itheta,
+ double[] directionError,
+ double[] shapeParameters,
+ int type,
+ int pid) {
+
+ addHits(hits);
+
+ // This will override the energy calculated from the hits, by design!
+ this.energy = energy;
+ this.energyError = energyError;
+
+ this.position = position;
+ this.positionError = positionError;
+
+ this.iphi = iphi;
+ this.itheta = itheta;
+
+ this.directionError = directionError;
+ this.shapeParameters = shapeParameters;
+
+ this.type = type;
+ this.pid = pid;
+
+ this.needsPropertyCalculation = false;
+ }
+
+ /**
+ * Basic constructor that takes a list of hits.
+ * It will apply the default energy calculation.
+ * @param hits The list of CalorimeterHits.
+ */
+ public BaseCluster(List<CalorimeterHit> hits) {
+ addHits(hits);
+ }
+
+ /**
+ *****************************************************
+ * Implementation of get methods from the interface. *
+ *****************************************************
+ */
+
+ /**
+ * Get the list of CalorimeterHits of this cluster.
+ * The hits are not necessarily unique in the list.
+ * @return The hits comprising the cluster.
+ */
+ public List<CalorimeterHit> getCalorimeterHits() {
+ return hits;
+ }
+
+ /**
+ * Get the clusters that are part of this cluster.
+ * @return The clusters comprising the cluster.
+ */
+ public List<Cluster> getClusters() {
+ return clusters;
+ }
+
+ /**
+ * Get the energy of the cluster, which by default will be the
+ * sum of the CalorimeterHit corrected energy values.
+ * @return The energy of the cluster.
+ */
+ public double getEnergy() {
+ return energy;
+ }
+
+ /**
+ * Get the energy error.
+ * @return The energy error.
+ */
+ public double getEnergyError() {
+ return energyError;
+ }
+
+ /**
+ * Get the individual hit contribution energies.
+ * This should be an array of the same size as the hit list.
+ * By default this array contains the hit's corrected energies,
+ * but the contributions may be set to different values
+ * on a per hit basis using the {@link #addHit(CalorimeterHit, double)}
+ * method.
+ * @return The individual hit contribution energies.
+ */
+ public double[] getHitContributions() {
+ double[] arrayCopy = new double[hitContributions.size()];
+ for (int i = 0; i < hitContributions.size(); i++) {
+ arrayCopy[i] = hitContributions.get(i);
+ }
+ return arrayCopy;
+ }
+
+ /**
+ * Get the list of subdetector energy contributions.
+ * The ordering and meaning of this array is unspecified by this class.
+ * @return The list of subdetector energy contributions.
+ */
+ public double[] getSubdetectorEnergies() {
+ return subdetectorEnergies;
+ }
+
+ /**
+ * Get a value defining the type of this cluster.
+ * @return The type of this cluster.
+ */
+ public int getType() {
+ return type;
+ }
+
+ /**
+ * Get the number of hits in the cluster, including hits in sub-clusters.
+ * Hits belonging to more than one cluster are counted once.
+ * @return The size of the cluster.
+ */
+ public int getSize() {
+ Set<CalorimeterHit> hitSet = new HashSet<CalorimeterHit>(this.hits);
+ for (Cluster cluster : clusters) {
+ hitSet.addAll(cluster.getCalorimeterHits());
+ }
+ int size = hitSet.size();
+ hitSet.clear();
+ return size;
+ }
+
+ /**
+ * Get the intrinsic phi direction of the cluster.
+ * @return The intrinsic phi direction of the cluster.
+ */
+ public double getIPhi() {
+ checkCalculateProperties();
+ return iphi;
+ }
+
+ /**
+ * Get the intrinsic theta direction of the cluster.
+ * @return The intrinsic theta direction of the cluster.
+ */
+ public double getITheta() {
+ checkCalculateProperties();
+ return itheta;
+ }
+
+ /**
+ * Get the direction error of the cluster as a double array of size 6.
+ * @return The direction error of the cluster.
+ */
+ public double[] getDirectionError() {
+ checkCalculateProperties();
+ return directionError;
+ }
+
+ /**
+ * Get the position of the cluster as a double array of size 3.
+ * @return The position of the cluster.
+ */
+ public double[] getPosition() {
+ checkCalculateProperties();
+ return position;
+ }
+
+ /**
+ * Get the position error of the cluster as a double array of size 6.
+ * @return The position error of the cluster.
+ */
+ public double[] getPositionError() {
+ checkCalculateProperties();
+ return positionError;
+ }
+
+ /**
+ * Get the shape parameters of the cluster as a double array of unspecified size.
+ * @return The shape parameters of the cluster.
+ */
+ public double[] getShape() {
+ checkCalculateProperties();
+ return shapeParameters;
+ }
+
+ /**
+ * Get the PDG ID of the particle hypothesis.
+ * @return The PID.
+ */
+ public int getParticleId() {
+ return pid;
+ }
+
+ /**
+ **********************************
+ * Implementation of set methods. *
+ **********************************
+ */
+
+ /**
+ * Set the position of the cluster.
+ * @param position The position of the cluster.
+ * @throws IllegalArgumentException if array is wrong size.
+ */
+ public void setPosition(double[] position) {
+ if (position.length != 3) {
+ throw new IllegalArgumentException("The position array argument has the wrong length: " + position.length);
+ }
+ this.position = position;
+ }
+
+ /**
+ * Set the position error of the cluster.
+ * @param positionError The position error of the cluster.
+ * @throws IllegalArgumentException if array is wrong size.
+ */
+ public void setPositionError(double[] positionError) {
+ if (positionError.length != 6) {
+ throw new IllegalArgumentException("The positionError array argument has the wrong length: " + position.length);
+ }
+ this.positionError = positionError;
+ }
+
+ /**
+ * Set the intrinsic phi of the cluster.
+ * @param iphi The intrinsic phi of the cluster.
+ */
+ public void setIPhi(double iphi) {
+ this.iphi = iphi;
+ }
+
+ /**
+ * Set the intrinsic theta of the cluster.
+ * @param iphi The intrinsic theta of the cluster.
+ */
+ public void setITheta(double itheta) {
+ this.itheta = itheta;
+ }
+
+ /**
+ * Set the direction error of the cluster.
+ * @param directionError The direction error of the cluster.
+ * @throws IllegalArgumentException if array is wrong size.
+ */
+ public void setDirectionError(double[] directionError) {
+ if (directionError.length != 6) {
+ throw new IllegalArgumentException("The directionError array argument has the wrong length: " + position.length);
+ }
+ this.directionError = directionError;
+ }
+
+ /**
+ * Set the shape parameters of the cluster.
+ * @param shapeParameters The shape parameters.
+ */
+ public void setShapeParameters(double[] shapeParameters) {
+ this.shapeParameters = shapeParameters;
+ }
+
+ /**
+ * Set the type of the cluster.
+ * @param type The type of the cluster.
+ */
+ public void setType(int type) {
+ this.type = type;
+ }
+
+ /**
+ * Get the PDG ID of the particle hypothesis.
+ * @return The PID.
+ */
+ public void setParticleId(int pid) {
+ this.pid = pid;
+ }
+
+ /**
+ * Set a total energy of this cluster, overriding any energy
+ * value that may have been automatically calculated from
+ * hit energies.
+ * @param energy The total energy of this cluster.
+ */
+ public void setEnergy(double energy) {
+ this.energy = energy;
+ }
+
+ /**
+ * Set the error on the energy measurement.
+ * @param energyError The error on the energ measurement.
+ */
+ public void setEnergyError(double energyError) {
+ this.energyError = energyError;
+ }
+
+ /**
+ * Set the subdetector energies.
+ * @param subdetectorEnergies The subdetector energies.
+ */
+ public void setSubdetectorEnergies(double[] subdetectorEnergies) {
+ this.subdetectorEnergies = subdetectorEnergies;
+ }
+
+ /**
+ ****************************************************
+ * Convenience methods for adding hits and clusters *
+ ****************************************************
+ */
+
+ /**
+ * Add a list of hits to the cluster.
+ * @param The list of hits to add.
+ */
+ public void addHits(List<CalorimeterHit> hits) {
+ for (CalorimeterHit hit : hits) {
+ addHit(hit);
+ }
+ }
+
+ /**
+ * Add a list of sub-clusters to the cluster.
+ * @param The list of clusters to add.
+ */
+ public void addClusters(List<Cluster> clusters) {
+ for (Cluster cluster : clusters) {
+ addCluster(cluster);
+ }
+ }
+
+ /**
+ * Add a sub-cluster to the cluster.
+ * @param cluster The cluster to add.
+ */
+ public void addCluster(Cluster cluster) {
+ clusters.add(cluster);
+ List<CalorimeterHit> clusterHits = cluster.getCalorimeterHits();
+ for (int i = 0; i < clusterHits.size(); i++) {
+ hits.add(clusterHits.get(i));
+ hitContributions.add(clusterHits.get(i).getCorrectedEnergy());
+ }
+ energy += cluster.getEnergy();
needsPropertyCalculation = true;
- }
-
- /**
- * Add a CalorimeterHit to the cluster from a cluster
- */
- public void addHitFromCluster(CalorimeterHit hit)
- {
- hits.add(hit);
- double hce = hit.getCorrectedEnergy();
- hit_energies.add(hce);
+ }
+
+ /**
+ * Add a hit to the cluster with default energy contribution.
+ * @param hit The hit to add.
+ */
+ public void addHit(CalorimeterHit hit) {
+ addHit(hit, hit.getCorrectedEnergy());
needsPropertyCalculation = true;
}
-
- /**
- * Remove a CalorimeterHit from the cluster
- */
- public void removeHit(CalorimeterHit hit)
- {
- int indx = hits.indexOf(hit);
+
+ /**
+ * Add a hit to the cluster with specified energy contribution.
+ * @param hit The hit to add.
+ * @param contribution The energy contribution of the hit [GeV].
+ */
+ public void addHit(CalorimeterHit hit, double contribution) {
+ hits.add(hit);
+ hitContributions.add(contribution);
+ energy += contribution;
+ needsPropertyCalculation = true;
+ }
+
+ /**
+ * Remove a hit from the cluster.
+ * @param hit The hit to remove.
+ */
+ public void removeHit(CalorimeterHit hit) {
+ int index = hits.indexOf(hit);
hits.remove(hit);
- double hre = hit.getRawEnergy();
-// raw_energy -= hre;
- double hce = hit.getCorrectedEnergy();
- IDDecoder idc = hit.getIDDecoder();
- idc.setID(hit.getCellID());
- int detector_index = idc.getSystemID();
- if((detector_index > 9 )|(detector_index < 0))detector_index = 0;
- corrected_energy -= hce;
- hit_energies.remove(indx);
- subdetector_raw_energies[detector_index] -= hre;
- subdetector_corrected_energies[detector_index] -= hce;
+ double hitEnergy = hit.getCorrectedEnergy();
+ energy -= hitEnergy;
+ hitContributions.remove(index);
needsPropertyCalculation = true;
}
- /**
- * Add a Cluster to the cluster
- */
- public void addCluster(Cluster cluster)
- {
- clusters.add(cluster);
- List<CalorimeterHit> chits = cluster.getCalorimeterHits();
- if(chits.size() > 0)
- {
- for(int i=0;i<chits.size();i++)
- {
- addHitFromCluster(chits.get(i));
- }
- needsPropertyCalculation = true;
- }
- corrected_energy += cluster.getEnergy();
- if(cluster instanceof BaseCluster)
- {
- BaseCluster bcl = (BaseCluster) cluster;
-// raw_energy += bcl.getRawEnergy();
- }
- double[] sde = cluster.getSubdetectorEnergies();
- for(int i=0;i<sde.length;i++)
- {
- subdetector_corrected_energies[i] += sde[i];
- subdetector_raw_energies[i] += sde[i];
- }
- }
-
- /**
- * Do an external calculation of Energy and set it
- */
- public void setEnergy(double Energy)
- {
- corrected_energy = Energy;
- }
-
- /**
- * Return the hits comprising the cluster as per interface
- */
- public List<CalorimeterHit> getCalorimeterHits()
- {
- return hits;
- }
-
- /**
- * Return the clusters comprising the cluster as per interface
- */
- public List<Cluster> getClusters()
- {
- return clusters;
- }
-
- /**
- * Return the corrected Energy of the cluster as per interface
- */
- public double getEnergy()
- {
- return corrected_energy;
- }
-
- public double getEnergyError()
- {
- return energyError;
- }
-
- public void setEnergyError(double energyError)
- {
- this.energyError = energyError;
- }
-
- /**
- * Return the corrected energy contribution of each hit
- * comprising the cluster as per interface
- */
- public double[] getHitContributions()
- {
- double[] earray = new double[hit_energies.size()];
- for(int i=0;i<hit_energies.size();i++)
- {
- earray[i] = hit_energies.get(i);
- }
- return earray;
- }
-
- /**
- * Return the phi direction of the cluster as per interface
- */
- public double getIPhi()
- {
- if(needsPropertyCalculation)
- {
+
+ /**
+ **************************************
+ * ClusterPropertyCalculator methods. *
+ **************************************
+ */
+
+ /**
+ * Calculate properties if needs property calculation.
+ */
+ void checkCalculateProperties() {
+ if (needsPropertyCalculation()) {
calculateProperties();
}
- return iphi;
- }
-
- /**
- * Return the theta direction of the cluster as per interface
- */
- public double getITheta()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return itheta;
- }
-
- /**
- * Return the direction error of the cluster as per interface
- */
- public double[] getDirectionError()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return directionError;
- }
-
- /**
- * Return the position of the cluster as per interface
- */
- public double[] getPosition()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return position;
- }
-
- /**
- * Return the position error of the cluster as per interface
- */
- public double[] getPositionError()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return positionError;
- }
-
- /**
- * Return the shape parameters of the cluster as per interface
- */
- public double[] getShape()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return shapeParameters;
- }
-
- /**
- * Return the subdetector energy contributions
- * of the cluster as per interface
- */
- public double[] getSubdetectorEnergies()
- {
- return subdetector_corrected_energies;
- }
-
- /**
- * Return a bit mask defining the type of cluster as per interface
- */
- public int getType()
- {
- return 0;
- }
-
-// /**
-// * Return the sum of the raw energies from the hits in the cluster
-// */
-// public double getRawEnergy()
-// {
-// return raw_energy;
-// }
-
- /**
- * Allow for recalculation of cluster properties, or
- * prevent recalculation after adding to cluster
- */
- public void setNeedsPropertyCalculation(boolean tf)
- {
- needsPropertyCalculation = tf;
- }
-
- /**
- * Override the default ClusterPropertyCalculator
- */
- public void setPropertyCalculator(ClusterPropertyCalculator cpc)
- {
- cluster_property_calculator = cpc;
- needsPropertyCalculation = true;
- }
-
- /**
- * Calculate the cluster properties from the hits
- */
- public void calculateProperties()
- {
- cluster_property_calculator.calculateProperties(hits);
- iphi = cluster_property_calculator.getIPhi();
- itheta = cluster_property_calculator.getITheta();
-
- // important: store info locally, not a reference to info in calculator
- double[] tmp = cluster_property_calculator.getPosition();
- for(int i=0; i<3; ++i) position[i] = tmp[i];
-
- double[] posErr = cluster_property_calculator.getPositionError();
- double[] dirErr = cluster_property_calculator.getDirectionError();
- for(int i=0; i<6; ++i) {
- positionError[i] = posErr[i];
- directionError[i] = dirErr[i];
- }
-
- // In general, assume any number of parameters needed
- double[] shapePar = cluster_property_calculator.getShapeParameters();
- shapeParameters = new double[ shapePar.length ];
- for(int i=0; i<shapePar.length; ++i) {
- shapeParameters[i] = shapePar[i];
- }
-
- needsPropertyCalculation = false;
- }
-
- /** Return the hits comprising the cluster, including hits in
- * subclusters, as per interface. Hits belonging to more than one
- * cluster/subcluster should be counted only once.
- */
- public int getSize()
- {
- Set<CalorimeterHit> hitSet = new HashSet<CalorimeterHit>(this.hits);
- for( Cluster clus : clusters )
- {
- hitSet.addAll( clus.getCalorimeterHits() );
- }
-
- // cleanup and return
- int size = hitSet.size();
- hitSet.clear();
- return size;
- }
-
- /**
- * String representation of this object.
- * @return A String representation of this object
- */
- public String toString()
- {
- double[] pos = this.getPosition();
- double rho = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
- double theta = Math.atan2( rho, pos[2]);
- double phi = Math.atan2( pos[1], pos[0] );
- if(phi<0) phi += 2*Math.PI;
- return "BaseCluster w/ "+this.getSize()+" hits at theta="+theta*180/Math.PI+", phi="+phi*180/Math.PI;
- }
-
- /**
- * Create and return a set of MCParticles from the hits of the cluster.
- * Only looks at SimCalorimeterHits, as CalorimeterHits don't have MC information.
- * @return The set of unique MCParticles associated with this cluster.
- */
- public Set<MCParticle> findUniqueMCParticlesFromHits()
- {
- Set<MCParticle> mcps = new HashSet<MCParticle>();
- for (CalorimeterHit hit : getCalorimeterHits())
- {
- if (hit instanceof SimCalorimeterHit)
- {
- SimCalorimeterHit simhit = (SimCalorimeterHit)hit;
- for (int i=0; i<simhit.getMCParticleCount(); i++)
- {
- mcps.add(simhit.getMCParticle(i));
- }
- }
- }
- return mcps;
+ }
+
+ /**
+ * Set a property calculator for computing position, etc.
+ * @param calc The property calculator.
+ */
+ public void setPropertyCalculator(ClusterPropertyCalculator calc) {
+ this.calc = calc;
+ }
+
+ /**
+ * True if property calculator is set.
+ * @return True if property calculator is set.
+ */
+ public boolean hasPropertyCalculator() {
+ return calc != null;
+ }
+
+ /**
+ * True if cluster is flagged as needed a property calculation.
+ * @return True if cluster needs property calculation.
+ */
+ public boolean needsPropertyCalculation() {
+ return needsPropertyCalculation;
+ }
+
+ /**
+ * Manually set whether the cluster needs property calculation.
+ * @param needsPropertyCalculation Whether the cluster needs property calculation.
+ */
+ public void setNeedsPropertyCalculation(boolean needsPropertyCalculation) {
+ this.needsPropertyCalculation = needsPropertyCalculation;
+ }
+
+ /**
+ *
+ * Calculate the properties of this cluster using the current
+ * <code>ClusterPropertyCalculator</code>. The calculated properties will be
+ * set on the following class variables:<br/>
+ * {@link #position}, {@link #positionError},
+ * {@link #iphi}, {@link #itheta}, {@link #directionError},
+ * and {@link #shapeParameters}. The cluster will then be flagged
+ * as not needing a property calculation in its current state.
+ */
+ public void calculateProperties() {
+ if (!hasPropertyCalculator()) {
+ throw new RuntimeException("No ClusterPropertyCalculator is set on this object.");
+ }
+ calc.calculateProperties(this);
+ setPosition(calc.getPosition());
+ setPositionError(calc.getPositionError());
+ setIPhi(calc.getIPhi());
+ setITheta(calc.getITheta());
+ setDirectionError(calc.getDirectionError());
+ setShapeParameters(calc.getShapeParameters());
+ setNeedsPropertyCalculation(false);
}
}
Modified: projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/ClusterPropertyCalculator.java
=============================================================================
--- projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/ClusterPropertyCalculator.java (original)
+++ projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/ClusterPropertyCalculator.java Wed Jan 7 17:03:40 2015
@@ -1,8 +1,7 @@
package org.lcsim.event.base;
-import java.util.ArrayList;
-import java.util.Collections;
import java.util.List;
+
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.Cluster;
@@ -11,41 +10,45 @@
*
* @author cassell
*/
-public interface ClusterPropertyCalculator
-{
+public interface ClusterPropertyCalculator {
- /**
- * Calculate properties from a CalorimeterHit list
- */
- public void calculateProperties(List<CalorimeterHit> hits);
+ /**
+ * Calculate properties from a CalorimeterHit list.
+ */
+ public void calculateProperties(List<CalorimeterHit> hits);
+
+ /**
+ * Calculate properties from a cluster.
+ */
+ public void calculateProperties(Cluster cluster);
- /**
- * Return position
- */
- public double[] getPosition();
-
- /**
- * Return position error
- */
- public double[] getPositionError();
-
- /**
- * Return phi direction
- */
- public double getIPhi();
-
- /**
- * Return theta direction
- */
- public double getITheta();
-
- /**
- * Return direction error
- */
- public double[] getDirectionError();
-
- /**
- * Return shape parameters
- */
- public double[] getShapeParameters();
+ /**
+ * Return position
+ */
+ public double[] getPosition();
+
+ /**
+ * Return position error
+ */
+ public double[] getPositionError();
+
+ /**
+ * Return phi direction
+ */
+ public double getIPhi();
+
+ /**
+ * Return theta direction
+ */
+ public double getITheta();
+
+ /**
+ * Return direction error
+ */
+ public double[] getDirectionError();
+
+ /**
+ * Return shape parameters
+ */
+ public double[] getShapeParameters();
}
Modified: projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/TensorClusterPropertyCalculator.java
=============================================================================
--- projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/TensorClusterPropertyCalculator.java (original)
+++ projects/lcsim/trunk/event-model/src/main/java/org/lcsim/event/base/TensorClusterPropertyCalculator.java Wed Jan 7 17:03:40 2015
@@ -1,26 +1,19 @@
package org.lcsim.event.base;
-import java.util.ArrayList;
-import java.util.Collections;
import java.util.List;
+
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.Cluster;
-// import org.lcsim.geometry.CalorimeterIDDecoder;
+
/**
- * Implementation of ClusterPropertyCalculator Interface using
- * inertia Tensor calculation to derive parameters
- *
+ * 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;
+public class TensorClusterPropertyCalculator extends AbstractClusterPropertyCalculator {
+
protected double[] mm_NE;
protected double[] mm_CE;
protected double[][] mm_PA;
@@ -28,26 +21,30 @@
/**
* 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)
- {
+ public TensorClusterPropertyCalculator() {
+ }
+
+ protected void reset() {
+ super.reset();
mm_NE = new double[3];
mm_CE = new double[3];
mm_PA = new double[3][3];
- for(int i=0;i<3;++i)
- {
+ }
+
+ public void calculateProperties(Cluster cluster) {
+ calculateProperties(cluster.getCalorimeterHits());
+ }
+
+ public void calculateProperties(List<CalorimeterHit> hits) {
+
+ reset();
+
+ 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.;}
+ for (int j = 0; j < 3; ++j) {
+ mm_PA[i][j] = 0.;
+ }
}
double Etot = 0.0;
double Exx = 0.0;
@@ -70,53 +67,51 @@
double M = 0.0;
double Det = 0.0;
int nhits = hits.size();
- for(int i=0;i<hits.size();i++)
- {
+ for (int i = 0; i < hits.size(); i++) {
CalorimeterHit hit = hits.get(i);
- // CalorimeterIDDecoder decoder = hit.getDecoder();
- // decoder.setID(hit.getCellID());
- // double[] pos = new double[3];
- // pos[0] = decoder.getX();
- // pos[1] = decoder.getY();
- // pos[2] = decoder.getZ();
+ // CalorimeterIDDecoder decoder = hit.getDecoder();
+ // decoder.setID(hit.getCellID());
+ // double[] pos = new double[3];
+ // pos[0] = decoder.getX();
+ // pos[1] = decoder.getY();
+ // pos[2] = decoder.getZ();
double[] pos = hit.getPosition();
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;
+ 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;
+ 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;
+ 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 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
@@ -124,180 +119,163 @@
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");
+ 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.);
}
- 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)
- {
+ if (mE1 < mE2) {
+ if (mE1 < mE3) {
E1 = mE1;
- if(mE2 < mE3)
- {
+ if (mE2 < mE3) {
E2 = mE2;
E3 = mE3;
- }
- else
- {
+ } else {
E2 = mE3;
E3 = mE2;
}
- }
- else
- {
+ } else {
E1 = mE3;
E2 = mE1;
E3 = mE2;
}
- }
- else
- {
- if(mE2 < mE3)
- {
+ } else {
+ if (mE2 < mE3) {
E1 = mE2;
- if(mE1 < mE3)
- {
+ if (mE1 < mE3) {
E2 = mE1;
E3 = mE3;
- }
- else
- {
+ } else {
E2 = mE3;
E3 = mE1;
}
- }
- else
- {
+ } else {
E1 = mE3;
E2 = mE2;
E3 = mE1;
}
}
- NE1 = E1/Etot;
- NE2 = E2/Etot;
- NE3 = E3/Etot;
+ 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;
- }
- }
+ // 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;
@@ -305,66 +283,57 @@
mm_CE[0] = CEx;
mm_CE[1] = CEy;
mm_CE[2] = CEz;
- for(int i=0;i<6;i++)
- {
+ for (int i = 0; i < 6; i++) {
positionError[i] = 0.;
directionError[i] = 0.;
shapeParameters[i] = 0.;
}
- for(int i=0;i<3;i++)
- {
+ 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]) ) ;
+ 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
- {
+ 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()
- {
+
+ public double[] getPosition() {
return position;
}
- public double[] getPositionError()
- {
+
+ public double[] getPositionError() {
return positionError;
}
- public double getIPhi()
- {
+
+ public double getIPhi() {
return iphi;
}
- public double getITheta()
- {
+
+ public double getITheta() {
return itheta;
}
- public double[] getDirectionError()
- {
+
+ public double[] getDirectionError() {
return directionError;
}
- public double[] getShapeParameters()
- {
+
+ public double[] getShapeParameters() {
return shapeParameters;
}
- public double[] getNormalizedEigenvalues()
- {
+
+ public double[] getNormalizedEigenvalues() {
return mm_NE;
}
- public double[][] getPrincipleAxis()
- {
+
+ public double[][] getPrincipleAxis() {
return mm_PA;
}
########################################################################
Use REPLY-ALL to reply to list
To unsubscribe from the LCDET-SVN list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCDET-SVN&A=1
