lcsim/src/org/lcsim/recon/pfa/identifier

diff -u -r1.2 -r1.3
--- MIPChargedParticleMaker.java	17 Aug 2006 20:55:37 -0000	1.2
+++ MIPChargedParticleMaker.java	3 Nov 2006 19:09:42 -0000	1.3
@@ -11,6 +11,11 @@

 import org.lcsim.mc.fast.tracking.ReconTrack;
 import org.lcsim.event.base.BaseReconstructedParticle;
 import org.lcsim.recon.ztracking.cheater.CheatTrack;

+import org.lcsim.event.MCParticle;
+import org.lcsim.geometry.Detector;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.recon.cluster.util.ClusterEnergyCalculator;
+

 
 /**
  * Given a list of MIP clusters and a list of tracks,

@@ -27,11 +32,10 @@

  * then the entire cluster is added to the ReconstructedParticle instead.
  * The parent must be unique.
  *

- * @version $Id: MIPChargedParticleMaker.java,v 1.2 2006/08/17 20:55:37 mcharles Exp $

+ * @version $Id: MIPChargedParticleMaker.java,v 1.3 2006/11/03 19:09:42 mcharles Exp $

  */
 

-public class MIPChargedParticleMaker
- extends Driver

+public class MIPChargedParticleMaker extends Driver

 {
     /** Simple constructor. */
     public MIPChargedParticleMaker() {

@@ -44,15 +48,35 @@

     public void setOutputMIPList(String name){ m_outputMIPListName = name; }
     public void setOutputParticleList(String name){ m_outputParticleListName = name; }
     public void setTrackMatcher(TrackClusterMatcher matcher) { m_matcher = matcher; }

-   

+
+    protected boolean m_checkEoverP = false;
+    /**
+     * Enable/disable a check on E/p. If enabled, a calibration must
+     * also be specified with <code>setCalibration()</code>.
+     */
+    public void setCheckEoverP(boolean check) { m_checkEoverP = check; }
+    /** 
+     * Specify an energy calibration. This is only used for the E/p
+     * check, if enabled. Note that setting a calibration doesn't
+     * automatically enable the E/p check.
+     */
+    public void setCalibration(ClusterEnergyCalculator calib) { m_calib = calib; }
+    protected ClusterEnergyCalculator m_calib = null;
+
+

     public void addClusterList(String inputName, String outputName) {
 	m_clusterLists.put(inputName, outputName);
     }
 

+    boolean m_debug = false;
+    public void setDebug(boolean debug) { m_debug = debug; }
+
+    protected EventHeader m_event = null;

     public void process(EventHeader event)
     {
 	super.process(event);

-	

+	m_event = event;
+

 	// Inputs:
 	List<Track> inputTrackList = event.get(Track.class, m_inputTrackListName);
 	List<Cluster> inputMIPList = event.get(Cluster.class, m_inputMIPListName);

@@ -148,7 +172,43 @@

 		matchedClusters.put(clus, part);
 	    }
 	}

-	outputParticleList.addAll(matchedClusters.values());

+
+	if (m_checkEoverP) {
+	    for (LocalReconstructedParticle part : matchedClusters.values()) {
+		boolean energyOK = checkEoverP(part);
+		if (energyOK) {
+		    outputParticleList.add(part);
+		} else {
+		    // No match => undo association
+		    // Remove track matching...
+		    List<Track> partTracks = part.getTracks();
+		    unmatchedTracks.addAll(partTracks);
+		    // Remove MIP matching...
+		    List<Cluster> partClusters = part.getClusters();
+		    for (Cluster partCluster : partClusters) {
+			List<Cluster> daughters = recursivelyFindSubClusters(partCluster);
+			for (Cluster dauClus : daughters) {
+			    if (inputMIPList.contains(dauClus)) {
+				unmatchedMIPs.add(dauClus);
+			    }
+			}
+		    }	  
+		    // Remove macro cluster matching...
+		    for (String inputListName : m_clusterLists.keySet()) {
+			String outputListName = m_clusterLists.get(inputListName);
+			List<Cluster> inputList = inputClusterLists.get(inputListName);
+			List<Cluster> outputList = outputClusterLists.get(outputListName);
+			for (Cluster partCluster : partClusters) {
+			    if (inputList.contains(partCluster) && !(outputList.contains(partCluster))) {
+				outputList.add(partCluster);
+			    }
+			}
+		    }
+		}
+	    }
+	} else {
+	    outputParticleList.addAll(matchedClusters.values());
+	}

 
 	// Write out
 	event.put(m_outputTrackListName, unmatchedTracks);

@@ -159,6 +219,11 @@

 	    List<Cluster> outputList = outputClusterLists.get(outputName);
 	    event.put(outputName, outputList);
 	}

+
+	if (m_debug) {
+	    System.out.println("MIPChargedParticleMaker: Read in "+inputTrackList.size()+" tracks and "+inputMIPList.size()+" MIP clusters; wrote out "
+			       +outputParticleList.size()+" matched particles, "+unmatchedTracks.size()+" unmatched tracks, "+unmatchedMIPs.size()+" unmatched MIP clusters");
+	}

     }
 
     /**

@@ -219,4 +284,46 @@

 	    this.setEnergy(energy);
 	}
     }

+
+    protected boolean checkEoverP(LocalReconstructedParticle part) {
+        // We don't expect an exact match due to resolution, energy lost
+        // to fragments etc., but a good portion of the energy should be
+        // in the cluster
+     
+        // Check energy and uncertainty from calorimeter:
+	List<Cluster> partClusters = part.getClusters();
+	double estimatedClusterEnergy = 0.0;
+	for (Cluster clus : partClusters) {
+          estimatedClusterEnergy += estimateClusterEnergy(clus);
+	}
+        double estimatedClusterEnergyUncertainty = 0.7 * Math.sqrt(estimatedClusterEnergy); // 70%/sqrt(E) for hadrons
+        
+        // Check energy from tracks
+	double estimatedTrackEnergy = 0.0;
+	for (Track tr : part.getTracks()) {
+	    double[] trackMomentum = tr.getMomentum();
+	    double trackMomentumMagSq = trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2];
+	    double trackMomentumMag = Math.sqrt(trackMomentumMagSq);
+	    double massPion   = 0.14;
+	    double energyReleaseIfPion = Math.sqrt(trackMomentumMagSq + massPion*massPion);
+	    estimatedTrackEnergy += energyReleaseIfPion;
+	}
+
+        double allowedVariation = 3.0; // 3sigma
+        boolean energyDiffOK_pion       = Math.abs((estimatedTrackEnergy - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty) < allowedVariation;
+        boolean energyDiffOK = energyDiffOK_pion;
+
+	if (m_debug) {
+	    System.out.println("DEBUG [MIP]: Checked E/P; found cluster E="+estimatedClusterEnergy
+			       +" +- "+estimatedClusterEnergyUncertainty
+			       +" and track E="+estimatedTrackEnergy
+			       +" => "+energyDiffOK);
+	}
+
+        return energyDiffOK;
+    }
+
+    private double estimateClusterEnergy(Cluster clus) {
+	return m_calib.getEnergy(clus);
+    }

 }