lcsim/src/org/lcsim/recon/pfa/identifier
diff -N SimpleChargedParticleMaker.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SimpleChargedParticleMaker.java 5 Jul 2006 00:10:29 -0000 1.1
@@ -0,0 +1,405 @@
+package org.lcsim.recon.pfa.structural.identifier;
+
+import java.util.*;
+import hep.physics.vec.*;
+import hep.physics.particle.Particle;
+
+import org.lcsim.event.Track;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.TensorClusterPropertyCalculator;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.base.BaseReconstructedParticle;
+
+/**
+ * Given lists of clusters and tracks, make a list of charged
+ * ReconstructedParticles.
+ *
+ * Currently, PID is done by cheating.
+ *
+ * @version $Id: SimpleChargedParticleMaker.java,v 1.1 2006/07/05 00:10:29 mcharles Exp $
+ */
+
+public class SimpleChargedParticleMaker extends Driver
+{
+ /** Simple constructor. */
+ public SimpleChargedParticleMaker() {}
+
+ public void setInputTrackList(String name) { m_inputTrackListName = name; }
+ public void setOutputTrackList(String name){ m_outputTrackListName = name; }
+ public void setInputMIPList(String name){ m_inputMIPListName = name; }
+ public void setInputClusterList(String name){ m_inputClusterListName = name; }
+ public void setOutputParticleList(String name){ m_outputParticleListName = name; }
+
+ public void process(EventHeader event)
+ {
+ initGeometry(event);
+
+ // Inputs:
+ List<Track> inputTrackList = event.get(Track.class, m_inputTrackListName);
+ List<Cluster> inputMIPList = event.get(Cluster.class, m_inputMIPListName);
+ List<Cluster> inputClusterList = event.get(Cluster.class, m_inputClusterListName);
+
+ // Outputs:
+ Map<Track,Cluster> matchedTracks = new HashMap<Track,Cluster>();
+ List<Track> outputTrackList = new Vector<Track>(); // initially empty
+ List<ReconstructedParticle> outputParticleList = new Vector<ReconstructedParticle>();
+
+ // Get the field strength using Steve McGill's code:
+ Detector det = event.getDetector();
+ double[] zero = {0, 0, 0};
+ double[] fieldStrength = det.getFieldMap().getField(zero);
+ HelixSwimmer swimmer = new HelixSwimmer(fieldStrength[2]);
+
+ // For each track, we'll try to extrapolate it to the ECAL surface.
+ // If we do, we store the appropriate value of alpha (distance to swim)
+ Map<Track, Double> successfullyExtrapolatedTracks = new HashMap<Track,Double>();
+
+ for (Track tr : inputTrackList) {
+ // Make a swimmer:
+ swimmer.setTrack(tr);
+ double alpha = Double.NaN;
+ // // Try swimming to the barrel:
+ double alphaBarrel = swimToBarrel(swimmer);
+ boolean validBarrel = false;
+ // Try swimming to the endcap:
+ double alphaEndcap = swimToEndcap(swimmer);
+ boolean validEndcap = false;
+ // Fixme: Here we should check that the track really does go all the
+ // way to the ECAL instead of stopping/decaying/interacting earlier.
+ // This used to be done with the checkDecayPoint() method in
+ // contrib.uiowa.structural.SwimToECAL
+ if (isValidBarrelIntercept(swimmer, alphaBarrel)) {
+ alpha = alphaBarrel;
+ validBarrel = true;
+ } else if (isValidEndcapIntercept(swimmer, alphaEndcap)) {
+ alpha = alphaEndcap;
+ validEndcap = true;
+ }
+
+ if ( ! Double.isNaN(alpha) ) {
+ // Found something.
+ if ( validEndcap || validBarrel) {
+ if ( !(validEndcap && validBarrel) ) {
+ successfullyExtrapolatedTracks.put(tr, new Double(alpha));
+ } else {
+ throw new AssertionError("Confusing intercept! barrel="+validBarrel+", endcap="+validEndcap);
+ }
+ } else {
+ throw new AssertionError("alpha is not NaN, but not a valid barrel or endcap intercept");
+ }
+ }
+ }
+
+ // See: structural.MatchHelixToCluster
+ for (Track tr : successfullyExtrapolatedTracks.keySet()) {
+ Double alpha = successfullyExtrapolatedTracks.get(tr);
+ Cluster matchedCluster = null; // This will be from inputClusterList
+ // First try to match it to a MIP...
+ Cluster matchedMIP = findMatchedMIP(tr, swimmer, alpha.doubleValue(), inputMIPList);
+ if (matchedMIP != null) {
+ // Matching MIP. Look up its parent cluster:
+ for (Cluster parent : inputClusterList) {
+ List<Cluster> clustersInsideParent = recursivelyFindSubClusters(parent);
+ if (clustersInsideParent.contains(matchedMIP)) {
+ // Found it
+ matchedCluster = parent;
+ break;
+ }
+ }
+ } else {
+ // If that didn't work, try any cluster...
+ matchedCluster = findMatchedCluster(tr, swimmer, alpha.doubleValue(), inputClusterList);
+ }
+
+ if (matchedCluster != null) {
+ // We found a match
+ matchedTracks.put(tr, matchedCluster);
+ } else {
+ // No match -- put the track into the list of unused tracks
+ outputTrackList.add(tr);
+ }
+ }
+
+ // Now, we need to make output particles for the tracks.
+ // We have to watch out for the special cases where:
+ // >1 track is matched to a cluster
+ Map<Cluster,LocalReconstructedParticle> matchedClusters = new HashMap<Cluster,LocalReconstructedParticle> ();
+ for (Track tr : matchedTracks.keySet()) {
+ Cluster clus = matchedTracks.get(tr);
+ if (matchedClusters.keySet().contains(clus)) {
+ // Already used this cluster in a particle => just add the track to that particle
+ LocalReconstructedParticle part = matchedClusters.get(clus);
+ part.addTrack(tr);
+ double[] trackMomentum = tr.getMomentum();
+ double trackMomentumMagSq = (trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2]);
+ double trackMass = 0.140; // Treat the second particle as a pion
+ if (tr instanceof ReconTrack) {
+ Particle truthParticle = ((ReconTrack)(tr)).getMCParticle();
+ if (truthParticle != null) {
+ trackMass = truthParticle.getMass();
+ }
+ }
+ double trackEnergy = Math.sqrt(trackMomentumMagSq + trackMass*trackMass);
+ part.setEnergy(part.getEnergy() + trackEnergy);
+ } else {
+ // This cluster hasn't been used yet -- initialize its particle
+ LocalReconstructedParticle part = new LocalReconstructedParticle();
+ part.addTrack(tr);
+ part.addCluster(clus);
+ double[] trackMomentum = tr.getMomentum();
+ double trackMomentumMagSq = (trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2]);
+ double mass = 0.140;
+ if (tr instanceof ReconTrack) {
+ Particle truthParticle = ((ReconTrack)(tr)).getMCParticle();
+ if (truthParticle != null) {
+ mass = truthParticle.getMass();
+ }
+ }
+ double energy = Math.sqrt(trackMomentumMagSq + mass*mass);
+ part.setEnergy(energy);
+ matchedClusters.put(clus, part);
+ }
+ }
+
+ outputParticleList.addAll(matchedClusters.values());
+
+ // Write out
+ event.put(m_outputTrackListName, outputTrackList);
+ event.put(m_outputParticleListName, outputParticleList);
+ }
+
+ protected double swimToBarrel(HelixSwimmer swimmer) {
+ // Look for a hit in the first layer of the ECAL barrel
+ return swimmer.getDistanceToRadius(m_ECAL_barrel_r);
+ }
+ protected double swimToEndcap(HelixSwimmer swimmer) {
+ // Look for a hit in the first layer of the ECAL endcap
+ return swimmer.getDistanceToZ(m_ECAL_endcap_z);
+ }
+ protected boolean isValidBarrelIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have -m_ECAL_barrel_z <= z <= +m_ECAL_barrel_z
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double z = intercept.z();
+ boolean zInRange = (z >= m_ECAL_barrel_zmin && z <= m_ECAL_barrel_zmax);
+ return zInRange;
+ }
+ protected boolean isValidEndcapIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have m_ECAL_endcap_rmin <= r <= m_ECAL_endcap_rmax
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double r = Math.sqrt(intercept.x()*intercept.x() + intercept.y()*intercept.y());
+ boolean rInRange = (r >= m_ECAL_endcap_rmin && r <= m_ECAL_endcap_rmax);
+ return rInRange;
+ }
+ protected List<Cluster> recursivelyFindSubClusters(Cluster clus)
+ {
+ List<Cluster> output = new Vector<Cluster>();
+ for (Cluster dau : clus.getClusters()) {
+ output.addAll(recursivelyFindSubClusters(dau));
+ }
+ output.add(clus);
+ return output;
+ }
+
+ protected void initGeometry(EventHeader event)
+ {
+ if (!m_init) {
+ Detector det = event.getDetector();
+ CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
+ CylindricalCalorimeter eme = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap"));
+ m_ECAL_barrel_zmin = emb.getZMin();
+ m_ECAL_barrel_zmax = emb.getZMax();
+ m_ECAL_barrel_r = emb.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_z = eme.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_rmin = eme.getInnerRadius();
+ m_ECAL_endcap_rmax = eme.getOuterRadius();
+ m_init = true;
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_zmin="+m_ECAL_barrel_zmin);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_zmax="+m_ECAL_barrel_zmax);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_r="+m_ECAL_barrel_r);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_z="+m_ECAL_endcap_z);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_rmin="+m_ECAL_endcap_rmin);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_rmax="+m_ECAL_endcap_rmax);
+ }
+ }
+
+ protected Cluster findMatchedMIP(Track tr, HelixSwimmer swimmer, double alpha, List<Cluster> mips)
+ {
+ // Find the track intercept and direction
+ swimmer.setTrack(tr);
+ Hep3Vector trackPoint = swimmer.getPointAtDistance(alpha);
+ // Obtain the unit vector giving the tangent:
+ double delta = 0.1;
+ if (alpha < 0) { delta *= -1.0; }
+ Hep3Vector aLittleFurther = swimmer.getPointAtDistance(alpha+delta);
+ Hep3Vector tangent = VecOp.unit(VecOp.sub(aLittleFurther, trackPoint));
+
+ List<Cluster> nearestMIPs = findNearestClusters(trackPoint, mips);
+ for (Cluster nearbyMIP : nearestMIPs) {
+ // Obtain geometrical info:
+ CalorimeterHit nearestHit = findNearestHit(trackPoint, nearbyMIP);
+ double separation = proximity(trackPoint, nearestHit);
+ int firstLayerHit = getLayer(nearestHit);
+ double unitDotProduct = findUnitDotProduct(tangent, nearbyMIP);
+ org.lcsim.geometry.Subdetector subdet = nearestHit.getSubdetector();
+ // Make cuts:
+ boolean goodSubDet = (subdet.getName().compareTo("EMBarrel")==0) || (subdet.getName().compareTo("EMEndcap")==0);
+ boolean goodFirstLayer = (firstLayerHit < 5);
+ boolean goodDotProduct = (Math.abs(unitDotProduct) > 0.85);
+ boolean goodSeparation = (separation < 50.0);
+ boolean foundMatch = goodSubDet && goodFirstLayer && goodDotProduct && goodSeparation;
+ if (foundMatch) {
+ // OK, made a good match
+ return nearbyMIP;
+ }
+ }
+ // No match
+ return null;
+ }
+
+ protected Cluster findMatchedCluster(Track tr, HelixSwimmer swimmer, double alpha, List<Cluster> clusters)
+ {
+ // Find the track intercept and direction
+ swimmer.setTrack(tr);
+ Hep3Vector trackPoint = swimmer.getPointAtDistance(alpha);
+
+ List<Cluster> nearestClusters = findNearestClusters(trackPoint, clusters);
+ for (Cluster nearbyCluster : nearestClusters) {
+ // Obtain geometrical info:
+ CalorimeterHit nearestHit = findNearestHit(trackPoint, nearbyCluster);
+ double separation = proximity(trackPoint, nearestHit);
+ int firstLayerHit = getLayer(nearestHit);
+ org.lcsim.geometry.Subdetector subdet = nearestHit.getSubdetector();
+ // Make cuts:
+ boolean goodSubDet = (subdet.getName().compareTo("EMBarrel")==0) || (subdet.getName().compareTo("EMEndcap")==0);
+ boolean goodFirstLayer = (firstLayerHit < 5);
+ boolean goodSeparation = (separation < 30.0);
+ boolean foundMatch = goodSubDet && goodFirstLayer && goodSeparation;
+ if (foundMatch) {
+ // OK, made a good match
+ return nearbyCluster;
+ }
+ }
+ // No match
+ return null;
+ }
+
+ protected List<Cluster> findNearestClusters(Hep3Vector point, List<Cluster> clusterList)
+ {
+ Map<Cluster,Double> mapClusterToDistance = new HashMap<Cluster, Double>();
+ List<Cluster> sortedListOfClusters = new Vector<Cluster>();
+ for (Cluster clus : clusterList) {
+ double dist = proximity(point, clus);
+ mapClusterToDistance.put(clus, new Double(dist));
+ sortedListOfClusters.add(clus);
+ }
+ Comparator<Cluster> comp = new CompareMapping<Cluster>(mapClusterToDistance);
+ Collections.sort(sortedListOfClusters, comp);
+ return sortedListOfClusters;
+ }
+ protected CalorimeterHit findNearestHit(Hep3Vector point, Cluster clus)
+ {
+ CalorimeterHit nearest = null;
+ double minDist = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ if (distance<minDist || nearest==null) {
+ nearest = hit;
+ }
+ }
+ return nearest;
+ }
+
+ protected double proximity(Hep3Vector point, Cluster clus) {
+ CalorimeterHit nearestHit = findNearestHit(point, clus);
+ return proximity(point, nearestHit);
+ }
+ protected double proximity(Hep3Vector point, CalorimeterHit hit) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ return distance;
+ }
+
+ protected double findUnitDotProduct(Hep3Vector tangent, Cluster clus)
+ {
+ // Find the cluster direction
+ BasicCluster copy = new BasicCluster();
+ copy.addCluster(clus);
+ TensorClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
+ copy.setPropertyCalculator(calc);
+ copy.calculateProperties();
+ double[][]axes = calc.getPrincipleAxis();
+ Hep3Vector clusterDir = new BasicHep3Vector(axes[0][0], axes[0][1], axes[0][2]);
+ // Get the dot product:
+ double unitDotProduct = VecOp.dot(tangent, clusterDir) / (tangent.magnitude() * clusterDir.magnitude());
+ return unitDotProduct;
+ }
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+
+ private class CompareMapping<T> implements Comparator<T> {
+ public CompareMapping(Map<T,Double> map) {
+ m_map = map;
+ }
+ public int compare(Object o1, Object o2) {
+ Cluster c1 = (Cluster) o1;
+ Cluster c2 = (Cluster) o2;
+ Double D1 = m_map.get(c1);
+ Double D2 = m_map.get(c2);
+ if (D1.equals(D2)) {
+ // Equal
+ return 0;
+ } else if (D1.doubleValue() < D2.doubleValue()) {
+ return -1;
+ } else {
+ return +1;
+ }
+ }
+ Map<T,Double> m_map;
+ }
+
+ protected boolean m_init = false;
+ protected double m_ECAL_barrel_zmin;
+ protected double m_ECAL_barrel_zmax;
+ protected double m_ECAL_barrel_r;
+ protected double m_ECAL_endcap_z;
+ protected double m_ECAL_endcap_rmin;
+ protected double m_ECAL_endcap_rmax;
+
+ String m_inputTrackListName;
+ String m_outputTrackListName;
+ String m_inputMIPListName;
+ String m_inputClusterListName;
+ String m_outputParticleListName;
+
+ // This should be replaced by a central implementation
+ private class LocalReconstructedParticle extends BaseReconstructedParticle
+ {
+ public void setEnergy(double e) {
+ _fourVec.setT(e);
+ }
+ public void setMomentum(Hep3Vector p3) {
+ _fourVec.setV3(_fourVec.t(), p3);
+ }
+ public void setReferencePoint(Hep3Vector p3) {
+ _referencePoint = p3;
+ }
+ public void setCharge(double q) {
+ _charge = q;
+ }
+ }
+}
+
+
lcsim/src/org/lcsim/recon/pfa/identifier
diff -N SimpleNeutralParticleMaker.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SimpleNeutralParticleMaker.java 5 Jul 2006 00:10:29 -0000 1.1
@@ -0,0 +1,139 @@
+package org.lcsim.recon.pfa.structural.identifier;
+
+import java.util.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.event.CalorimeterHit;
+
+/**
+ * Given lists of clusters, make a list of neutral ReconstructedParticles.
+ *
+ * Currently, PID is predetermined and the energy calibration is ad-hoc.
+ *
+ * @version $Id: SimpleNeutralParticleMaker.java,v 1.1 2006/07/05 00:10:29 mcharles Exp $
+ */
+
+public class SimpleNeutralParticleMaker extends Driver
+{
+ /**
+ * Constructor.
+ *
+ * @param pdg The neutral particles are assumed to be of this type when calculating the 4-momentum.
+ */
+ public SimpleNeutralParticleMaker(int pdg) {
+ setParticleID(pdg);
+ }
+
+ // Configure
+ public void setInputClusterList(String name){ m_inputClusterListName = name; }
+ public void setOutputParticleList(String name){ m_outputParticleListName = name; }
+ public void setParticleID(int pdg) {
+ hep.physics.particle.properties.ParticlePropertyProvider mgr =
+ hep.physics.particle.properties.ParticlePropertyManager.getParticlePropertyProvider();
+ hep.physics.particle.properties.ParticleType type = mgr.get(pdg);
+ double m_mass = type.getMass();
+ }
+
+ // Process one event
+ public void process(EventHeader event)
+ {
+ // Input, output:
+ List<Cluster> inputClusterList = event.get(Cluster.class, m_inputClusterListName);
+ List<ReconstructedParticle> outputParticleList = new Vector<ReconstructedParticle>();
+
+ for (Cluster clus : inputClusterList) {
+ LocalReconstructedParticle part = new LocalReconstructedParticle();
+ part.addCluster(clus);
+ double energy = estimateClusterEnergy(clus); // clumsy, not a real calibration...
+ part.setEnergy(energy);
+ // Set the other particle properties that are needed to render
+ // properly in the event display.
+ part.setMomentum(computeMomentum(energy, clus));
+ part.setReferencePoint(new BasicHep3Vector(0,0,0));
+ part.setCharge(0);
+ // Add to the output list
+ outputParticleList.add(part);
+ }
+
+ event.put(m_outputParticleListName, outputParticleList);
+ }
+
+ private Hep3Vector computeMomentum(double energy, Cluster clus) {
+ // Where is the cluster?
+ Hep3Vector pos = new BasicHep3Vector(clus.getPosition());
+ Hep3Vector unitDirection = VecOp.unit(pos);
+ // Now, what's the momentum?
+ // p^2 = E^2 - m^2
+ double momentumSquared = energy*energy - m_mass*m_mass;
+ double momentumMagnitude = Math.sqrt(momentumSquared); // may complain if -ve...
+ Hep3Vector momentum = VecOp.mult(momentumMagnitude, unitDirection);
+ return momentum;
+ }
+
+ String m_inputClusterListName;
+ String m_outputParticleListName;
+ double m_mass;
+
+ // This should be replaced by a central implementation
+ private class LocalReconstructedParticle extends BaseReconstructedParticle
+ {
+ public void setEnergy(double e) {
+ _fourVec.setT(e);
+ }
+ public void setMomentum(Hep3Vector p3) {
+ _fourVec.setV3(_fourVec.t(), p3);
+ }
+ public void setReferencePoint(Hep3Vector p3) {
+ _referencePoint = p3;
+ }
+ public void setCharge(double q) {
+ _charge = q;
+ }
+ }
+
+ // Ad-hoc calibration, not general
+ protected double estimateClusterEnergy(Cluster clus) {
+ 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
+ clusterEnergySum += hit.getCorrectedEnergy();
+ //clusterEnergySum += hit.getRawEnergy() * 86.8;
+ } else if (name.compareTo("EMEndcap") == 0) {
+ // EM endcap -- OK
+ clusterEnergySum += hit.getCorrectedEnergy();
+ //clusterEnergySum += hit.getRawEnergy() * 86.8;
+ } else if (name.compareTo("HADBarrel") == 0) {
+ // Had barrel -- count
+ 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
+ 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;
+ }
+
+}