4 added files
lcsim/src/org/lcsim/contrib/uiowa/util
diff -N EnergyCalibration.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EnergyCalibration.java 16 Dec 2005 21:13:33 -0000 1.1
@@ -0,0 +1,23 @@
+package util;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.MCParticle;
+
+/**
+ * An interface for energy calibrations.
+ */
+
+public interface EnergyCalibration {
+
+ /**
+ * Return the corrected energy for this cluster.
+ */
+ public double energy(Cluster clus);
+
+ /**
+ * Return the corrected energy contributed to this cluster by this particle.
+ */
+ public double energyContributedByParticle(Cluster clus, MCParticle part);
+
+
+}
lcsim/src/org/lcsim/contrib/uiowa/util
diff -N HitCountAssociator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HitCountAssociator.java 16 Dec 2005 21:13:33 -0000 1.1
@@ -0,0 +1,92 @@
+package util;
+
+import java.util.List;
+import java.util.Comparator;
+import java.util.Collections;
+import java.util.Set;
+import java.util.HashSet;
+import java.util.Map;
+import java.util.HashMap;
+
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.EventHeader;
+
+public class HitCountAssociator extends EnergyAssociator
+{
+ public HitCountAssociator(EventHeader event) { super(event) ; }
+
+ /**
+ * Given a Cluster, returns a list of all MCParticles that
+ * contributed, sorted by amount of hits contributed.
+ */
+ public List<MCParticle> associateClusterToMCParticles(Cluster clus)
+ {
+ Set<MCParticle> contributingParticles = new HashSet<MCParticle>();
+ List<CalorimeterHit> hits = clus.getCalorimeterHits();
+ for (CalorimeterHit hit : hits) {
+ SimCalorimeterHit simHit = (SimCalorimeterHit) hit;
+ int nContributingParticles = simHit.getMCParticleCount();
+ for (int i=0; i<nContributingParticles; i++) {
+ contributingParticles.add(simHit.getMCParticle(i));
+ }
+ }
+ LocalSortedMap<MCParticle> map = new LocalSortedMap<MCParticle>();
+ for (MCParticle part : contributingParticles) {
+ map.put(part, new Double(getHitsContributedToCluster(part, clus)));
+ }
+ List<MCParticle> outputList = map.getSortedList();
+ // Crosscheck
+ MCParticle previous = null;
+ for (MCParticle tmpPart : outputList) {
+ if (previous != null) {
+ int prevCount = getHitsContributedToCluster(previous, clus);
+ int thisCount = getHitsContributedToCluster(tmpPart, clus);
+ if (thisCount > prevCount) {
+ String bling = new String();
+ bling += "BUG in sorting! prevCount="+prevCount+" and thisCount="+thisCount+", but map.get(prev)=";
+ bling += map.get(previous);
+ bling += " and map.get(this)=";
+ bling += map.get(tmpPart);
+ throw new AssertionError(bling);
+ }
+ }
+ previous = tmpPart;
+ }
+ return outputList;
+ }
+
+ /**
+ * Given an MCParticle, returns a list of all Clusters to which it
+ * contributes, sorted by amount of energy contributed.
+ */
+ public List<Cluster> associateMCParticleToClusters(MCParticle part)
+ {
+ // Map from cluster to amount of enery contributed
+ LocalSortedMap<Cluster> map = new LocalSortedMap<Cluster>();
+ for (Cluster clus : m_clusters) {
+ int hits = getHitsContributedToCluster(part, clus);
+ if (hits > 0) {
+ map.put(clus, new Double(hits));
+ }
+ }
+ return map.getSortedList();
+ }
+
+ public int getHitsContributedToCluster(MCParticle part, Cluster clus)
+ {
+ int hitCount = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ SimCalorimeterHit simHit = (SimCalorimeterHit) hit;
+ int nContributingParticles = simHit.getMCParticleCount();
+ for (int i=0; i<nContributingParticles; i++) {
+ if (part == simHit.getMCParticle(i)) {
+ hitCount++;
+ }
+ }
+ }
+ return hitCount;
+ }
+}
lcsim/src/org/lcsim/contrib/uiowa/util
diff -N OldEnergyCalibration.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ OldEnergyCalibration.java 16 Dec 2005 21:13:33 -0000 1.1
@@ -0,0 +1,88 @@
+package util;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.CalorimeterHit;
+
+/**
+ * Mat's ad-hoc energy calibration for sidaug05
+ */
+
+public class OldEnergyCalibration implements EnergyCalibration {
+
+ /**
+ * Constructor.
+ */
+ public OldEnergyCalibration() {
+ //
+ }
+
+ /**
+ * Return the corrected energy for this cluster.
+ */
+ public double energy(Cluster clus) {
+ return energy(clus, null, true, true);
+ }
+
+
+ /**
+ * Take just the hits in this cluster deposited by this particle and compute
+ * the corrected energy for them.
+ */
+ public double energyContributedByParticle(Cluster clus, MCParticle part) {
+ return energy(clus, part, true, true);
+ }
+
+ public double energyECAL(Cluster clus) {
+ return energy(clus, null, true, false);
+ }
+
+ public double energyHCAL(Cluster clus) {
+ return energy(clus, null, false, true);
+ }
+
+ protected double energy(Cluster clus, MCParticle part, boolean useECAL, boolean useHCAL)
+ {
+ double clusterEnergySum = 0.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) {
+ // EM barrel -- OK
+ if (useECAL) {
+ clusterEnergySum += hit.getCorrectedEnergy();
+ }
+ } else if (name.compareTo("EMEndcap") == 0) {
+ // EM endcap -- OK
+ if (useECAL) {
+ clusterEnergySum += hit.getCorrectedEnergy();
+ }
+ } else if (name.compareTo("HADBarrel") == 0) {
+ // Had barrel -- count
+ if (useHCAL) {
+ double rawEnergy = hit.getRawEnergy();
+ boolean goodHit = (hit.getTime() < 100);
+ double energyPerHit = 0.115;
+ if (goodHit) {
+ clusterEnergySum += energyPerHit;
+ }
+ }
+ } else if (name.compareTo("HADEndcap") == 0) {
+ // Had endcap -- count
+ if (useHCAL) {
+ double rawEnergy = hit.getRawEnergy();
+ boolean goodHit = (hit.getTime() < 100);
+ double energyPerHit = 0.115;
+ if (goodHit) {
+ clusterEnergySum += energyPerHit;
+ }
+ }
+ } else {
+ throw new AssertionError("DEBUG: Found non-calorimeterhit in calorimeter '"+subdet+"' with name '"+subdet.getName()+"'");
+ }
+ }
+
+ return clusterEnergySum;
+ }
+}
lcsim/src/org/lcsim/contrib/uiowa/util
diff -N RonEnergyCalibration.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RonEnergyCalibration.java 16 Dec 2005 21:13:33 -0000 1.1
@@ -0,0 +1,280 @@
+package util;
+
+import java.util.List;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.geometry.CylindricalSubdetector;
+import org.lcsim.geometry.Detector;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+
+/**
+ * Ron's energy calibrations.
+ *
+ * Currently, only the sidaug05 is in there.
+ */
+
+public class RonEnergyCalibration implements EnergyCalibration {
+
+ /**
+ * Constructor. EventHeader is required to obtain calibration.
+ */
+ public RonEnergyCalibration(EventHeader event) {
+ m_event = event;
+ }
+
+ /**
+ * Return the corrected energy for this cluster.
+ */
+ public double energy(Cluster clus)
+ {
+ Detector det = m_event.getDetector();
+ String detName = det.getName();
+ if (detName.compareTo("sidaug05") == 0) {
+ return energy_sidaug05(clus);
+ } else {
+ throw new AssertionError("Detector '"+detName+"' not handled yet.");
+ }
+ }
+
+ /**
+ * Take just the hits in this cluster deposited by this particle and compute
+ * the corrected energy for them. NB! This is non-linear. For example if
+ * a cluster contains equal contributions from particles A and B then
+ * energyContributedByParticle_sidaug05(clus,A)=energyContributedByParticle_sidaug05(clus,B)
+ * but
+ * energyContributedByParticle_sidaug05(clus,A)+energyContributedByParticle_sidaug05(clus,B)>energyContributedByParticle_sidaug05(clus)
+ */
+ public double energyContributedByParticle(Cluster clus, MCParticle part)
+ {
+ if (part==null) { throw new AssertionError("Error: MCParticle 'part' is null"); }
+
+ Detector det = m_event.getDetector();
+ String detName = det.getName();
+ if (detName.compareTo("sidaug05") == 0) {
+ return energyContributedByParticle_sidaug05(clus, part);
+ } else {
+ throw new AssertionError("Detector '"+detName+"' not handled yet.");
+ }
+ }
+
+ /**
+ * Corrected energy for cluster (with sidaug05 calibration)
+ */
+ protected double energy_sidaug05(Cluster clus)
+ {
+ return energyContributedByParticle_sidaug05(clus, null);
+ }
+
+
+ /**
+ * Energy contribution from this particle, using the sidaug05 calibration.
+ * Warning: non-linear.
+ *
+ * If part==null, count all hits (irrespective of which particle they came from)
+ */
+ protected double energyContributedByParticle_sidaug05(Cluster clus, MCParticle part)
+ {
+ if (isHadronicCheating(clus)) {
+ // Use Ron's neutral hadron calibration:
+ checkHits(clus, part);
+ double EMEraw = m_rawEnergyECAL;
+ int nEMhits = m_rawHitsECAL;
+ int nHADhits = m_rawHitsHCAL;
+ double sangle = getSinPolarAngle(clus);
+ return energy_sidaug05(EMEraw, nEMhits, nHADhits, sangle);
+ } else {
+ // Use default E/M calibration:
+ double defaultCalibrationEnergySum = 0.0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ defaultCalibrationEnergySum += hit.getCorrectedEnergy();
+ }
+ return defaultCalibrationEnergySum;
+ }
+ }
+
+ /**
+ * Does this hit have a contribution from the given MCParticle?
+ */
+ 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;
+ }
+
+ /**
+ * How much raw energy does the MC Particle contribute to the hit?
+ */
+ 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;
+ }
+
+ /**
+ * Internal calibration for neutral hadrons
+ */
+ protected double energy_sidaug05(double EMEraw, int nEMhits, int nHADhits, double sangle)
+ {
+ // Here is Ron's snippet:
+ 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;
+ return (Emeas + offset);
+ }
+
+ /**
+ * Does this look like a hadronic shower?
+ * If so, use the hadron calibration.
+ * If not, use an E/M calibration.
+ */
+ protected boolean isHadronicCheating(Cluster clus)
+ {
+ HitCountAssociator assoc = new HitCountAssociator(m_event);
+ List<MCParticle> assocParticles = assoc.associateClusterToMCParticles(clus);
+ MCParticle dominant = assocParticles.iterator().next();
+ int pdg = dominant.getPDGID();
+ if (pdg==11 || pdg==-11 || pdg==22) {
+ // Clearly electromagnetic
+ return false;
+ } else {
+ // Hadronic (or other non-E/M)
+ return true;
+ }
+ }
+
+ /**
+ * Look at the hits and figure out what was deposited by this particle.
+ * If part==null, then count all hits.
+ */
+ protected void checkHits(Cluster clus, MCParticle part)
+ {
+ m_rawEnergyECAL = 0.0;
+ m_rawHitsECAL = 0;
+ m_rawHitsHCAL = 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
+ m_rawHitsECAL++;
+ if (part == null) {
+ m_rawEnergyECAL += hit.getRawEnergy();
+ } else {
+ SimCalorimeterHit simHit = (SimCalorimeterHit) hit;
+ m_rawEnergyECAL += hitMatchEnergy(part, simHit);
+ }
+ } else if (name.compareTo("HADBarrel") == 0 || name.compareTo("HADEndcap") == 0) {
+ // HCAL
+ if (hit.getTime() < 100) {
+ // Within first 100 ns => OK
+ m_rawHitsHCAL++;
+ }
+ } else {
+ throw new AssertionError("DEBUG: Found non-calorimeterhit in calorimeter '"+subdet+"' with name '"+subdet.getName()+"'");
+ }
+ }
+ }
+ }
+
+ /**
+ * Utility routine to get sin(polar angle)
+ */
+ protected double getSinPolarAngle(Cluster clus)
+ {
+ // Get the energy-weighted position:
+ 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; }
+ }
+ // Normalize:
+ for (int i=0; i<3; i++) {
+ energyWeightedPositionSum[i] /= eRawSum ;
+ }
+
+ // Get the sine of the polar angle:
+ double rSquared = energyWeightedPositionSum[0]*energyWeightedPositionSum[0] + energyWeightedPositionSum[1]*energyWeightedPositionSum[1];
+ double magSquared = rSquared + energyWeightedPositionSum[2]*energyWeightedPositionSum[2];
+ double sangle = Math.sqrt(rSquared / magSquared);
+ return sangle;
+ }
+
+ protected double m_rawEnergyECAL;
+ protected int m_rawHitsECAL;
+ protected int m_rawHitsHCAL;
+ protected EventHeader m_event;
+
+
+ // Utility classes.
+ // Not used right now, but I anticipate needing to count hits/energy from
+ // just certain places/particles/hits in a general way in the future.
+
+ /*
+ * private interface Cut {
+ * public boolean useThisHit(CalorimeterHit hit, MCParticle part);
+ * }
+ * private class CutRequireECAL {
+ * public CutRequireECAL() {}
+ * public boolean useThisHit(CalorimeterHit hit, MCParticle part) {
+ * org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ * if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ * String name = subdet.getName();
+ * return (name.compareTo("EMBarrel") == 0 || name.compareTo("EMEndcap") == 0);
+ * }
+ * }
+ * private class CutRequireHCAL {
+ * public CutRequireHCAL() {}
+ * public boolean useThisHit(CalorimeterHit hit, MCParticle part) {
+ * org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ * if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ * String name = subdet.getName();
+ * return (name.compareTo("HADBarrel") == 0 || name.compareTo("HADEndcap") == 0);
+ * }
+ * }
+ */
+}
CVSspam 0.2.8