lcsim/src/org/lcsim/contrib/uiowa

diff -u -r1.3 -r1.4
--- ReclusterDTreeDriver.java	18 Dec 2007 00:50:41 -0000	1.3
+++ ReclusterDTreeDriver.java	19 Dec 2007 22:44:47 -0000	1.4
@@ -386,7 +386,7 @@

 	}
 
 	List<Cluster> tmp = new Vector<Cluster>();

-	printStatus("FINAL STATUS:", tracksSortedByMomentum, allSharedClusters, newMapTrackToShowerComponents, newMapShowerComponentToTrack, newMapTrackToThreshold, newMapTrackToTolerance, tmp, tmp, linkableClusters, tmp, tmp, newMapTrackToVetoedAdditions);

+	printStatus("FINAL STATUS:", tracksSortedByMomentum, allSharedClusters, newMapTrackToShowerComponents, newMapShowerComponentToTrack, newMapTrackToThreshold, newMapTrackToTolerance, photons, mips, clumps, treesWithNoStructure, seedLeftoverHitClusters, newMapTrackToVetoedAdditions);

 
 	// Outputs
 

@@ -395,6 +395,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;

@@ -411,6 +412,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) {

@@ -437,6 +439,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 (...)

@@ -545,8 +579,14 @@

 	    outputNeutralHadronParticleList.add(part);
 	}
 

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

 	// Write out
 	m_event.put(m_outputParticleListName, outputParticleList);

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

 	m_event = null;
     }