lcsim/src/org/lcsim/contrib/uiowa

diff -u -r1.11 -r1.12
--- ReclusterDriver.java	17 Dec 2007 23:47:06 -0000	1.11
+++ ReclusterDriver.java	19 Dec 2007 22:37:58 -0000	1.12
@@ -30,7 +30,7 @@

   *
   * This version is PRELIMINARY.
   *

-  * @version $Id: ReclusterDriver.java,v 1.11 2007/12/17 23:47:06 mcharles Exp $

+  * @version $Id: ReclusterDriver.java,v 1.12 2007/12/19 22:37:58 mcharles Exp $

   * @author Mat Charles
   */
 

@@ -826,6 +826,7 @@

 	List<ReconstructedParticle> outputChargedParticleList = new Vector<ReconstructedParticle>();
 	List<ReconstructedParticle> outputPhotonParticleList = new Vector<ReconstructedParticle>();
 	List<ReconstructedParticle> outputNeutralHadronParticleList = new Vector<ReconstructedParticle>();

+	List<ReconstructedParticle> outputChargedParticleListWithEoverPveto = new Vector<ReconstructedParticle>();

 
 	// PID info
 	int pdg_piplus = 211;

@@ -842,6 +843,7 @@

 	ParticleType type_K0 = ParticlePropertyManager.getParticlePropertyProvider().get(pdg_K0);
 	BaseParticleID pid_K0 = new BaseParticleID(type_K0);
 	double mass_K0 = type_K0.getMass();

+	double windowEoverPveto = 2.5;

 
 	// Charged tracks
 	for (Track tr : tracksSortedByMomentum) {

@@ -849,7 +851,8 @@

 	    // write out charged shower
 	    BaseReconstructedParticle part = new BaseReconstructedParticle();
 	    for (Cluster clus : showerComponents) { part.addCluster(clus); }

-            // [FIXME: ADD FUZZY SIMCALORIMETERHITS]

+	    Cluster sharedHitClus = makeClusterOfSharedHits(showerComponents, allSharedClusters);
+	    if (sharedHitClus.getCalorimeterHits().size()>0) { part.addCluster(sharedHitClus); }

 	    part.addTrack(tr);
 	    part.setCharge(tr.getCharge());
 	    // Get the particle's three-momentum.

@@ -867,6 +870,38 @@

 	    part.setParticleIdUsed(pid_piplus);
 	    outputParticleList.add(part);
 	    outputChargedParticleList.add(part);

+	    // Now consider E/p veto
+	    double clusterEnergy = energy(showerComponents, allSharedClusters);
+	    double sigma = 0.7 * Math.sqrt(trackMomentumMag);
+	    if (trackMomentumMag < 1.0) { 
+		// Don't lowball uncertainty for low-momentum tracks
+		sigma = 0.7; 
+	    } 
+	    double normalizedResidual = (clusterEnergy-trackMomentumMag)/sigma;
+	    if (normalizedResidual > -windowEoverPveto && normalizedResidual < windowEoverPveto) {
+		// Passes E/p check
+		outputChargedParticleListWithEoverPveto.add(part);
+	    } else {
+		// Fails E/p check
+		BaseReconstructedParticle calPart = new BaseReconstructedParticle();
+		Cluster clus = makeCombinedCluster(showerComponents);
+		calPart.addCluster(clus);
+		if (sharedHitClus.getCalorimeterHits().size()>0) { calPart.addCluster(sharedHitClus); }
+		double clusterMomentumMagSq = clusterEnergy*clusterEnergy - mass_piplus*mass_piplus;
+		if (clusterMomentumMagSq<0.0) { clusterMomentumMagSq = 0.0; }
+		Hep3Vector calPos = new BasicHep3Vector(clus.getPosition());
+		Hep3Vector calThreeMomentum = VecOp.mult(Math.sqrt(clusterMomentumMagSq), VecOp.unit(calPos)); // assume it comes from the IP
+		HepLorentzVector calFourMomentum = new BasicHepLorentzVector(clusterEnergy, calThreeMomentum);
+		calPart.set4Vector(calFourMomentum);
+		calPart.setReferencePoint(0,0,0);
+		calPart.setCharge(0);
+		// Set the PID and mass
+		calPart.addParticleID(pid_K0);
+		calPart.setParticleIdUsed(pid_K0);
+		calPart.setMass(mass_piplus);
+		// Add to the output list
+		outputChargedParticleListWithEoverPveto.add(calPart);
+	    }

 	}
 
 	// Photons

@@ -876,7 +911,9 @@

 		// write out photon
 		BaseReconstructedParticle part = new BaseReconstructedParticle();
 		part.addCluster(clus);

-                // [FIXME: ADD FUZZY SIMCALORIMETERHITS]

+                // Add fuzzy hits
+		Cluster sharedHitClus = makeClusterOfSharedHits(clus, allSharedClusters);
+		if (sharedHitClus.getCalorimeterHits().size()>0) { part.addCluster(sharedHitClus); }

 		double clusterEnergy = energy(clus, allSharedClusters, m_photonCalib);
 		Hep3Vector pos = new BasicHep3Vector(clus.getPosition());
 		Hep3Vector threeMomentum = VecOp.mult(clusterEnergy, VecOp.unit(pos)); // assume it comes from the IP

@@ -972,7 +1009,9 @@

 	    // write out neutral hadron
 	    BaseReconstructedParticle part = new BaseReconstructedParticle();
 	    part.addCluster(clus);

-            // [FIXME: ADD FUZZY SIMCALORIMETERHITS]

+	    // Add fuzzy hits
+	    Cluster sharedHitClus = makeClusterOfSharedHits(clus, allSharedClusters);
+	    if (sharedHitClus.getCalorimeterHits().size()>0) { part.addCluster(sharedHitClus); }

 	    double clusterEnergy = energy(clus, allSharedClusters, m_neutralCalib);
 	    double clusterMomentumMagSq = clusterEnergy*clusterEnergy - mass_K0*mass_K0;
 	    if (clusterMomentumMagSq<0.0) { clusterMomentumMagSq = 0.0; }

@@ -1000,7 +1039,14 @@

 		System.out.println("DEBUG: Looking at neutral cluster with measured energy = "+clusterEnergy+". Dominant particle is "+truthPDG+" with p="+truthMomentum+". Cluster has effic="+effic+", purity="+purity);
 	    }
 	}

+
+	List<ReconstructedParticle> outputParticleListWithEoverPveto = new Vector<ReconstructedParticle>();
+	outputParticleListWithEoverPveto.addAll(outputPhotonParticleList);
+	outputParticleListWithEoverPveto.addAll(outputNeutralHadronParticleList);
+	outputParticleListWithEoverPveto.addAll(outputChargedParticleListWithEoverPveto);
+

 	m_event.put(m_outputParticleListName, outputParticleList);

+	m_event.put(m_outputParticleListName+"_withEoverPveto", outputParticleListWithEoverPveto);

 
 	makePlotsOfEoverP(outputChargedParticleList, newMapTrackToShowerComponents, allSharedClusters);