lcsim/src/org/lcsim/contrib/uiowa/structural/likelihood
diff -u -r1.4 -r1.5
--- MiscUtilities.java 18 Oct 2005 22:32:57 -0000 1.4
+++ MiscUtilities.java 16 Dec 2005 21:09:28 -0000 1.5
@@ -16,9 +16,16 @@
import org.lcsim.geometry.CylindricalSubdetector;
//import org.lcsim.geometry.Subdetector;
import org.lcsim.util.swim.Line;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.EventHeader;
-class MiscUtilities
+import util.HitCountAssociator;
+
+public class MiscUtilities
{
+ //MiscUtilities() {}
+
/**
* Get the minimum distance between hits in the two clusters.
* If one or both of the clusters is empty, will return NaN.
@@ -45,26 +52,34 @@
static protected double distance(Cluster clus, CalorimeterHit hit)
{
- // Loop over hits...
+ Hep3Vector point = new BasicHep3Vector(hit.getPosition());
+ return distance(clus, point);
+ }
+
+ static protected double distance(CalorimeterHit hit1, CalorimeterHit hit2)
+ {
+ Hep3Vector point2 = new BasicHep3Vector(hit2.getPosition());
+ return distance(hit1, point2);
+ }
+
+ static protected double distance(Cluster clus, Hep3Vector point)
+ {
boolean firstCheck = true;
double minDistance = Double.NaN; // Will stay NaN if clus is empty
List<CalorimeterHit> hits = clus.getCalorimeterHits();
for (CalorimeterHit hitInCluster : hits) {
- double dist = distance(hit, hitInCluster);
+ double dist = distance(hitInCluster, point);
if (firstCheck || dist<minDistance) {
minDistance = dist;
firstCheck = false;
}
}
-
return minDistance;
}
- static protected double distance(CalorimeterHit hit1, CalorimeterHit hit2)
- {
- Hep3Vector vect1 = new BasicHep3Vector(hit1.getPosition());
- Hep3Vector vect2 = new BasicHep3Vector(hit2.getPosition());
- Hep3Vector displacement = VecOp.sub(vect1, vect2);
+ static protected double distance(CalorimeterHit hit, Hep3Vector point) {
+ Hep3Vector hitPoint = new BasicHep3Vector(hit.getPosition());
+ Hep3Vector displacement = VecOp.sub(hitPoint, point);
return displacement.magnitude();
}
@@ -188,4 +203,213 @@
}
return matchingLayer;
}
+
+ // Blah... belongs in Line.
+ static protected double findDOCAToPoint(Hep3Vector linePoint, Hep3Vector lineDirection, Hep3Vector point)
+ {
+ // The first line is kind of ugly.
+ Hep3Vector displacement = new BasicHep3Vector(point.x() - linePoint.x(), point.y() - linePoint.y(), point.z() - linePoint.z());
+ double dotProduct = VecOp.dot(displacement, lineDirection);
+ double doca = Math.sqrt(displacement.magnitudeSquared() - dotProduct*dotProduct);
+ return doca;
+ }
+
+ /*
+ static protected boolean isHadronic(Cluster clus, EventHeader event)
+ {
+ throw new AssertionError("Not implemented, apparently");
+ //HitCountAssociator assoc = new HitCountAssociator(event);
+ //List<MCParticle> assocParticles = assoc.associateClusterToMCParticles(clus);
+ //MCParticle dominant = assocParticles.iterator().next();
+ //return ( Math.abs(dominant.getCharge()) < 0.001);
+ }
+
+ static protected double getCorrectedEnergy_sidaug05(Cluster clus, EventHeader event)
+ {
+ // Need to figure out the raw EM energy, the number of EM and HAD hits
+ double EMEraw = 0.0;
+ int nEMhits = 0;
+ int nHADhits = 0;
+ int nRejectedHitsHAD = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ String name = subdet.getName();
+ if (name.compareTo("EMBarrel") == 0 || name.compareTo("EMEndcap") == 0) {
+ // ECAL
+ nEMhits++;
+ EMEraw += hit.getRawEnergy();
+ } else if (name.compareTo("HADBarrel") == 0 || name.compareTo("HADEndcap") == 0) {
+ // HCAL
+ if (hit.getTime() < 100) {
+ // Within first 100 ns => OK
+ nHADhits++;
+ } else {
+ nRejectedHitsHAD++;
+ }
+ } else {
+ throw new AssertionError("DEBUG: Found non-calorimeterhit in calorimeter '"+subdet+"' with name '"+subdet.getName()+"'");
+ }
+ }
+
+ // Also need the "sangle", i.e. sin(theta)
+ double sangle = Double.NaN;
+ Hep3Vector[] posAndDir = MiscUtilities.getPositionAndDirection(clus);
+ Hep3Vector pos = null;
+ if (posAndDir != null) {
+ pos = posAndDir[0];
+ } else {
+ // Too few hits -- calculate energy-weighted position by hand.
+ double eRawSum = 0.0;
+ double[] energyWeightedPositionSum = new double[3];
+ for (int i=0; i<3; i++) { energyWeightedPositionSum[i] = 0.0; }
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ double eRaw = hit.getRawEnergy();
+ eRawSum += eRaw;
+ double[] hitPos = hit.getPosition();
+ for (int i=0; i<3; i++) { energyWeightedPositionSum[i] += hitPos[i] * eRaw; }
+ }
+ for (int i=0; i<3; i++) { energyWeightedPositionSum[i] /= eRawSum ; }
+ pos = new BasicHep3Vector(energyWeightedPositionSum);
+ }
+ double mag = pos.magnitude();
+ double r = Math.sqrt(pos.x()*pos.x() + pos.y()*pos.y());
+ sangle = r/mag;
+
+ // Now, is this an E/M shower or a hadronic shower?
+ boolean isHadronic = MiscUtilities.isHadronic(clus, event);
+
+ if (isHadronic) {
+ // Hadronic -- use Ron's routine
+ double correctedEnergy = MiscUtilities.getCorrectedEnergy_sidaug05(EMEraw, nEMhits, nHADhits, sangle);
+ return correctedEnergy;
+ } else {
+ // Electromagnetic -- use default calibration
+ double defaultCalibrationEnergySum = 0.0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ defaultCalibrationEnergySum += hit.getCorrectedEnergy();
+ }
+ return defaultCalibrationEnergySum;
+ }
+ }
+
+ static protected double getCorrectedEnergy_sidaug05(double EMEraw, int nEMhits, int nHADhits, double sangle)
+ {
+ final double HADsf = 8.81;
+ final double HADoffset = 4.;
+ final double EMlin = .033897;
+ final double EMquad = -.00000781;
+ final double EMelin = 81.418;
+ final double EMequad = -8.65;
+ final int EMncut = 200;
+ final double AngleCoef = 1.;
+ final double offset = 0.;
+
+ double Emeas = 0.;
+ double Emhad = 0.;
+ double Emem = 0.;
+ if(nEMhits < EMncut) {
+ Emem = nEMhits*EMlin + nEMhits*nEMhits*EMquad;
+ } else {
+ Emem = EMEraw*EMelin + EMEraw*EMEraw*EMequad;
+ }
+ double angcorr = 1. + AngleCoef*(Math.sqrt(sangle) -1.);
+ Emhad = (nHADhits+HADoffset)/HADsf/angcorr;
+ Emeas = Emem + Emhad;
+ //System.out.println("DEBUG: Emem="+Emem+", Emhad="+Emhad);
+ return (Emeas + offset);
+ }
+ */
+
+ static protected double correctedEnergyInClusterECAL(MCParticle part, Cluster clus)
+ {
+ double energyECAL = 0.0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean useThisHit = true;
+ if (part != null) {
+ useThisHit = hitMatch(part, (SimCalorimeterHit)(hit));
+ }
+ if (useThisHit) {
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ String name = subdet.getName();
+ if (name.compareTo("EMBarrel") == 0 || name.compareTo("EMEndcap") == 0) {
+ // ECAL
+ SimCalorimeterHit simHit = (SimCalorimeterHit) hit;
+ double rawEnergy = hitMatchEnergy(part, simHit);
+ double rawEnergyTotal = hit.getRawEnergy();
+ double frac = (rawEnergy / rawEnergyTotal);
+ energyECAL += frac * hit.getCorrectedEnergy();
+ }
+ }
+ }
+ return energyECAL;
+ }
+
+ static protected int countHitsInClusterHCAL(MCParticle part, Cluster clus)
+ {
+ int countHitsHCAL = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean useThisHit = true;
+ if (part != null) {
+ useThisHit = hitMatch(part, (SimCalorimeterHit)(hit));
+ }
+ if (useThisHit) {
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ String name = subdet.getName();
+ if (name.compareTo("HADBarrel") == 0 || name.compareTo("HADEndcap") == 0) {
+ // HCAL
+ if (hit.getTime() < 100) {
+ // Within first 100 ns => OK
+ countHitsHCAL++;
+ }
+ }
+ }
+ }
+ return countHitsHCAL;
+ }
+
+ static protected boolean hitMatch(MCParticle part, SimCalorimeterHit hit) {
+ int nContributingParticles = hit.getMCParticleCount();
+ for (int i=0; i<nContributingParticles; i++) {
+ MCParticle hitPart = hit.getMCParticle(i);
+ if (part == hitPart) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ static protected double hitMatchEnergy(MCParticle part, SimCalorimeterHit hit) {
+ int nContributingParticles = hit.getMCParticleCount();
+ for (int i=0; i<nContributingParticles; i++) {
+ MCParticle hitPart = hit.getMCParticle(i);
+ if (part == hitPart) {
+ return hit.getContributedEnergy(i);
+ }
+ }
+ return 0.0;
+ }
+
+ static protected int countHitsInCluster(MCParticle part, Cluster clus)
+ {
+ int countHits = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ if (hit instanceof SimCalorimeterHit) {
+ SimCalorimeterHit simHit = (SimCalorimeterHit) (hit);
+ int numContributingParticles = simHit.getMCParticleCount();
+ for (int i=0; i<numContributingParticles; i++) {
+ MCParticle contributingParticle = simHit.getMCParticle(i);
+ if (part == contributingParticle) {
+ // Match
+ countHits++;
+ }
+ }
+ } else {
+ throw new AssertionError("Non-simulated calorimeter hit");
+ }
+ }
+ return countHits;
+ }
}