lcsim/src/org/lcsim/contrib/uiowa

diff -u -r1.11 -r1.12
--- ReclusterDTreeDriver.java	11 Mar 2008 19:27:40 -0000	1.11
+++ ReclusterDTreeDriver.java	20 Mar 2008 02:29:18 -0000	1.12
@@ -40,6 +40,8 @@

     protected boolean m_useNewMipFinder = true;
     protected boolean m_useOldMipFinder = true;
     protected boolean m_clusterAsJets = true;

+    protected boolean m_ignorePunchThroughTracksForJets = true;
+    protected boolean m_useTracksThatDontReachCalorimeter = true;

 
     protected boolean m_useNewPhotonID = true;
 

@@ -576,6 +578,27 @@

 	unmatchedTracks.addAll(trackList);
 	unmatchedTracks.removeAll(tracksMatchedToClusters.keySet());
 	m_histo_unmatchedTracksMultiplicity.fill(unmatchedTracks.size());

+	List<Track> unmatchedTracksThatDontReachCalorimeter = new Vector<Track>();
+	for (Track tr : unmatchedTracks) {
+	    LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
+	    debugTrackMatch.process(m_event);
+	    Cluster debugMatchedCluster = debugTrackMatch.matchTrackToCluster(tr, allMatchableClusters);
+	    if (debugMatchedCluster != null) {
+		// This can happen sometimes if it pointed to a photon that was broken up severely.
+		// In any case, it clearly pointed to the calorimeter so we shouldn't
+		// add on the track momentum (that would be double-counting)
+	    } else {
+		Hep3Vector interceptPoint = debugTrackMatch.extendToLayer(0);
+		if (interceptPoint == null) {
+		    // No valid extrap to calorimeter
+		    unmatchedTracksThatDontReachCalorimeter.add(tr);
+		} else {
+		    // Reached calorimeter, but no match => we'll treat as NH or similar
+		}
+	    }
+	}
+
+	// debug

 	{
 	    double momentumSum = 0.0;
 	    for (Track tr : unmatchedTracks) {

@@ -960,13 +983,17 @@

 		double sigma = 0.7 * Math.sqrt(trackMomentum);
 		double normalizedResidual = (clusterEnergy-trackMomentum)/sigma;
 		boolean punchThrough = isPunchThrough(showerComponents, allSharedClusters);

-		if (normalizedResidual < -1.0 && !punchThrough) {

+		if (normalizedResidual < -1.0) {

 		    // FIXME: This cut-off shouldn't be hard-coded

-		    tracksWithEoverPTooLow.add(tr);

+		    if (!m_ignorePunchThroughTracksForJets || !punchThrough) {
+			tracksWithEoverPTooLow.add(tr);
+		    }

 		}

-		if (normalizedResidual < -3.0 && !punchThrough) {

+		if (normalizedResidual < -3.0) {

 		    // FIXME: This cut-off shouldn't be hard-coded

-		    tracksWithEoverPMuchTooLow.add(tr);

+		    if (!m_ignorePunchThroughTracksForJets || !punchThrough) {
+			tracksWithEoverPMuchTooLow.add(tr);
+		    }

 		}
 	    }
 

@@ -1873,32 +1900,33 @@

 	    }
 	    // write out neutral hadron
 	    BaseReconstructedParticle part = makeNeutralHadron(clus, allSharedClusters);

-	    /*
-	    BaseReconstructedParticle part = new BaseReconstructedParticle();
-	    part.addCluster(clus);
-	    // Add fuzzy hits
-	    Cluster sharedHitClus = makeClusterOfSharedHits(clus, allSharedClusters);
-	    if (sharedHitClus.getCalorimeterHits().size()>0) { part.addCluster(sharedHitClus); }
-	    double clusterEnergy = energy(clus, allSharedClusters, m_neutralCalib); // CHECK: Is this OK?
-	    double clusterMomentumMagSq = clusterEnergy*clusterEnergy - mass_K0*mass_K0;
-	    if (clusterMomentumMagSq<0.0) { clusterMomentumMagSq = 0.0; }
-	    Hep3Vector pos = new BasicHep3Vector(clus.getPosition());
-	    Hep3Vector threeMomentum = VecOp.mult(Math.sqrt(clusterMomentumMagSq), VecOp.unit(pos)); // assume it comes from the IP
-	    HepLorentzVector fourMomentum = new BasicHepLorentzVector(clusterEnergy, threeMomentum);
-	    part.set4Vector(fourMomentum);
-	    part.setReferencePoint(0,0,0);
-	    part.setCharge(0);
-	    // Set the PID and mass
-	    part.addParticleID(pid_K0);
-	    part.setParticleIdUsed(pid_K0);
-	    part.setMass(mass_K0);
-	    */

 	    // Add to the output list
 	    outputParticleList.add(part);
 	    outputNeutralHadronParticleList.add(part);
 	    if (m_debug) { debugPrintNeutralHadronParticleInfo(clus, allSharedClusters); }
 	}
 

+	// Tracks that didn't reach the ECAL
+	List<ReconstructedParticle> chargedParticlesThatDontReachCalorimeter = new Vector<ReconstructedParticle>();
+	for (Track tr : unmatchedTracksThatDontReachCalorimeter) {
+	    BaseReconstructedParticle part = new BaseReconstructedParticle();
+	    part.setCharge(tr.getCharge());
+	    Hep3Vector trackMomentum = momentum(tr);
+	    double trackMomentumMag = trackMomentum.magnitude();
+	    double currentParticleMass = mass_piplus;
+	    double trackEnergySq = trackMomentumMag*trackMomentumMag + currentParticleMass*currentParticleMass;
+	    HepLorentzVector fourMomentum = new BasicHepLorentzVector(Math.sqrt(trackEnergySq), trackMomentum);
+	    part.set4Vector(fourMomentum);
+	    part.setMass(mass_piplus);
+	    part.setReferencePoint(new BasicHep3Vector(tr.getReferencePoint()));
+	    part.addParticleID(pid_piplus);
+	    part.setParticleIdUsed(pid_piplus);
+	    chargedParticlesThatDontReachCalorimeter.add(part);
+	}
+	if (m_useTracksThatDontReachCalorimeter) {
+	    outputParticleList.addAll(chargedParticlesThatDontReachCalorimeter);
+	}
+

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

@@ -2280,6 +2308,7 @@

 	m_event.put(m_outputParticleListName+"_debugEoverP_oldCalib", outputChargedParticleListWithClusterEnergy_oldCalib);
 	m_event.put(m_outputParticleListName+"_debugEoverP_punchThrough_oldCalib", outputChargedParticleListWithClusterEnergy_punchThrough_oldCalib);
 	m_event.put(m_outputParticleListName+"_debugEoverP_noPunchThrough_oldCalib", outputChargedParticleListWithClusterEnergy_noPunchThrough_oldCalib);

+	m_event.put(m_outputParticleListName+"_dontReachCalorimeter", chargedParticlesThatDontReachCalorimeter);

 
 	m_event.put(m_outputParticleListName+"_photonLikePhotons", photonLikePhotons_particles);
 	m_event.put(m_outputParticleListName+"_chargedHadronLikePhotons", chargedHadronLikePhotons_particles);