lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.49 -r1.50
--- ReclusterDTreeDriver.java 9 Sep 2008 22:35:50 -0000 1.49
+++ ReclusterDTreeDriver.java 22 Sep 2008 23:59:25 -0000 1.50
@@ -35,7 +35,7 @@
* in this package, which uses the implementation in
* org.lcsim.recon.cluster.directedtree developed by NIU).
*
- * @version $Id: ReclusterDTreeDriver.java,v 1.49 2008/09/09 22:35:50 mcharles Exp $
+ * @version $Id: ReclusterDTreeDriver.java,v 1.50 2008/09/22 23:59:25 tjkim Exp $
* @author Mat Charles <[log in to unmask]>
*/
@@ -47,6 +47,8 @@
protected boolean m_cheatOnHadronsMisidentifiedAsPhotons = false;
protected boolean m_cheatOnPhotons = false;
+ protected boolean m_muonDebug = false;
+ protected boolean m_electronDebug = false;
protected boolean m_photonDebug = false;
protected boolean m_photonSplitDebug = false;
protected boolean m_jetDebug = false;
@@ -130,6 +132,7 @@
FindSubClusters clusDriverHCAL = new FindSubClusters("DTreeClustersHCAL", m_newMipFinderRadiusHCAL, m_minHitsToBeTreatedAsClusterHCAL, m_removePoorQualityMips, "OldMipsInsideTreesHCAL", "NewMipsInsideTreesHCAL", "ClumpsInsideTreesHCAL", "BlocksInsideTreesHCAL", "LeftoverHitsInsideTreesHCAL", "MapTreeToTargetsHCAL", "MapSharedToTreeHCAL");
FindSubClusters clusDriverMCAL = new FindSubClusters("DTreeClustersMCAL", m_newMipFinderRadiusMCAL, m_minHitsToBeTreatedAsClusterMCAL, m_removePoorQualityMips, "OldMipsInsideTreesMCAL", "NewMipsInsideTreesMCAL", "ClumpsInsideTreesMCAL", "BlocksInsideTreesMCAL", "LeftoverHitsInsideTreesMCAL", "MapTreeToTargetsMCAL", "MapSharedToTreeMCAL");
FindSubClusters clusDriverFCAL = new FindSubClusters("DTreeClustersFCAL", m_newMipFinderRadiusFCAL, m_minHitsToBeTreatedAsClusterFCAL, m_removePoorQualityMips, "OldMipsInsideTreesFCAL", "NewMipsInsideTreesFCAL", "ClumpsInsideTreesFCAL", "BlocksInsideTreesFCAL", "LeftoverHitsInsideTreesFCAL", "MapTreeToTargetsFCAL", "MapSharedToTreeFCAL");
+ FindSubClusters clusDriverMCALforMuonID = new FindSubClusters("DTreeClustersMCALforMuonID", m_newMipFinderRadiusMCAL, m_minHitsToBeTreatedAsClusterMCAL, m_removePoorQualityMips, "OldMipsInsideTreesMCALforMuonID", "NewMipsInsideTreesMCALforMuonID", "ClumpsInsideTreesMCALforMuonID", "BlocksInsideTreesMCALforMuonID", "LeftoverHitsInsideTreesMCALforMuonID", "MapTreeToTargetsMCALforMuonID", "MapSharedToTreeMCALforMuonID");
if (m_oldMipFinderCrossesTrees) {
clusDriverECAL.enableBarrelEndcapCrossing("EcalBarrDigiHits", "EcalBarrDigiHitsNearBoundary", "EcalEndcapDigiHits", "EcalEndcapDigiHitsNearBoundary");
clusDriverHCAL.enableBarrelEndcapCrossing("HcalBarrDigiHits", "HcalBarrDigiHitsNearBoundary", "HcalEndcapDigiHits", "HcalEndcapDigiHitsNearBoundary");
@@ -138,6 +141,7 @@
clusDriverECAL.setNNrange(1,1,1);
clusDriverHCAL.setNNrange(2,2,1);
clusDriverMCAL.setNNrange(2,2,1);
+ clusDriverMCALforMuonID.setNNrange(2,2,1);
clusDriverFCAL.setNNrange(1,1,1);
}
add(clusDriverECAL);
@@ -161,6 +165,12 @@
add(new TransientFlagDriver("ClumpsInsideTreesMCAL"));
add(new TransientFlagDriver("BlocksInsideTreesMCAL"));
add(new TransientFlagDriver("LeftoverHitsInsideTreesMCAL"));
+ add(clusDriverMCALforMuonID);
+ add(new TransientFlagDriver("OldMipsInsideTreesMCALforMuonID"));
+ add(new TransientFlagDriver("NewMipsInsideTreesMCALforMuonID"));
+ add(new TransientFlagDriver("ClumpsInsideTreesMCALforMuonID"));
+ add(new TransientFlagDriver("BlocksInsideTreesMCALforMuonID"));
+ add(new TransientFlagDriver("LeftoverHitsInsideTreesMCALforMuonID"));
}
if (m_useFcal) {
add(clusDriverFCAL);
@@ -270,6 +280,10 @@
mips.addAll(mipsOld);
mips.addAll(mipsNew);
+ //For muon ID
+ List<Cluster> mipsMuon = new Vector<Cluster>();
+ mipsMuon.addAll(event.get(Cluster.class, "OldMipsInsideTreesMCALforMuonID"));
+ mipsMuon.addAll(event.get(Cluster.class, "NewMipsInsideTreesMCALforMuonID"));
if (m_debug) {
System.out.println("Found "+mips.size()+" mips, "+clumps.size()+" clumps, "+photons.size()+" photons, "+leftoverHitClusters.size()+" leftover-hit-clusters, "+treesWithNoStructure.size()+" large DTrees with no structure, and "+trackList.size()+" tracks in event.");
}
@@ -327,6 +341,7 @@
List<Cluster> photonLikePhotons = new Vector<Cluster>();
List<Cluster> electronClusters = new Vector<Cluster>();
List<Track> electronTracks = new Vector<Track>();
+ List<Track> muonTracks = new Vector<Track>();
photonHandling(photons, electronClusters, chargedHadronLikePhotons, modifiedPhotonClusters, photonLikePhotons, trackList, electronTracks, clustersMatchedToTracks, tracksMatchedToClusters);
// Resume track matching
@@ -496,6 +511,14 @@
event.put("ShowerFinderMips", preShowerMips);
event.getMetaData(preShowerMips).setTransient(true);
+ //Muon finder
+ scanForMuons(trackList, MapTrkToMIP, mipsMuon, muonTracks);
+ event.put("MuonList", muonTracks);
+ if(m_muonDebug) { aida.cloud1D("muon/muon size",100000).fill(muonTracks.size());}
+ //Electron write out
+ event.put("ElectronList", electronTracks);
+ if(m_electronDebug) aida.cloud1D("electron size",100000).fill(electronTracks.size());
+
// Figure out whether tracks were uniquely matched or not:
Set<Track> uniquelyMatchedTracks = new HashSet<Track>();
Set<Track> ambiguouslyMatchedTracks = new HashSet<Track>();
@@ -1147,18 +1170,19 @@
}
ReassignClustersAlgorithm algorithm = null;
+ String mapName;
+ if (m_useSteveMipsForConeScoring) {
+ mapName = "TrackMipMap";
+ } else {
+ mapName = "ShowerFinderMapTrackToMip";
+ }
+ Map<Track, BasicCluster> mapTrackToMIP = (Map<Track, BasicCluster>) (m_event.get(mapName));
+ MIPGeometryHandler geomHandler = new LayerBasedMIPGeometryHandler(mapTrackToMIP, m_findCluster);
+ //MIPGeometryHandler geomHandler = new HelixTangentMIPGeometryHandler(mapTrackToMIP, m_findCluster);
+
if (m_useSimpleConeForReassignment) {
algorithm = new ConeReassignmentAlgorithm(1.00, m_findCluster);
} else {
- String mapName;
- if (m_useSteveMipsForConeScoring) {
- mapName = "TrackMipMap";
- } else {
- mapName = "ShowerFinderMapTrackToMip";
- }
- Map<Track, BasicCluster> mapTrackToMIP = (Map<Track, BasicCluster>) (m_event.get(mapName));
- MIPGeometryHandler geomHandler = new LayerBasedMIPGeometryHandler(mapTrackToMIP, m_findCluster);
- //MIPGeometryHandler geomHandler = new HelixTangentMIPGeometryHandler(mapTrackToMIP, findCluster);
algorithm = new ConeMIPReassignmentAlgorithm(geomHandler, 800.0, Math.PI*0.5);
}
if (m_fixSingleTracksWithCone) {
@@ -1167,13 +1191,13 @@
// Only process tracks that aren't part of a jet:
Set<Track> jetOfTrack = mapTrackToJet.get(tr);
if (jetOfTrack == null) {
- checkTrackForReassignments(tr, newMapTrackToShowerComponents, newMapShowerComponentToTrack, allSharedClusters, unassignedClusters, newMapTrackToTolerance.get(tr), algorithm, tweakedTracksMatchedToClusters);
+ checkTrackForReassignments(tr, newMapTrackToShowerComponents, newMapShowerComponentToTrack, allSharedClusters, unassignedClusters, newMapTrackToTolerance.get(tr), algorithm);
}
}
}
if (m_fixJetsWithCone) {
for (Set<Track> jet : jets) {
- checkJetForReassignments(jet, newMapJetToShowerComponents, newMapShowerComponentToJet, allSharedClusters, unassignedClusters, algorithm, newMapShowerComponentToTrack, tweakedTracksMatchedToClusters, newMapTrackToShowerComponents);
+ checkJetForReassignments(jet, newMapJetToShowerComponents, newMapShowerComponentToJet, allSharedClusters, unassignedClusters, algorithm, newMapTrackToShowerComponents);
}
}
@@ -1184,7 +1208,7 @@
// Outputs
makeParticlesAndWriteOut(trackList, tracksSortedByMomentum, unmatchedTracksThatDontReachCalorimeter, mapOrigTrackToTweakedTrack,
tracksMatchedToClusters, clustersMatchedToTracks,
- electronTracks, modifiedPhotonClusters,electronClusters, photonLikePhotons, chargedHadronLikePhotons,
+ electronTracks, muonTracks, modifiedPhotonClusters,electronClusters, photonLikePhotons, chargedHadronLikePhotons,
seedHadronLikePhotonClusters, nonSeedHadronLikePhotonClusters, nonSeedPhotonLikePhotonClusters,
newMapTrackToShowerComponents, newMapShowerComponentToTrack,
linkableClusters,
@@ -1232,11 +1256,10 @@
List<SharedClusterGroup> allSharedClusters,
Set<Cluster> unassignedClusters,
double toleranceOfTrack,
- ReassignClustersAlgorithm algorithm,
- Map<Track, Cluster> tweakedTracksMatchedToClusters)
+ ReassignClustersAlgorithm algorithm)
{
Set<Cluster> showerComponents = newMapTrackToShowerComponents.get(tr);
- List<Cluster> reassignedClusters = reassignClustersToTrack(tr, showerComponents, unassignedClusters, allSharedClusters, toleranceOfTrack, algorithm, newMapShowerComponentToTrack, tweakedTracksMatchedToClusters, newMapTrackToShowerComponents);
+ List<Cluster> reassignedClusters = reassignClustersToTrack(tr, showerComponents, unassignedClusters, allSharedClusters, toleranceOfTrack, algorithm, newMapTrackToShowerComponents);
if (reassignedClusters != null && reassignedClusters.size()>0) {
for (Cluster clus : reassignedClusters) {
showerComponents.add(clus);
@@ -1252,12 +1275,10 @@
List<SharedClusterGroup> allSharedClusters,
Set<Cluster> unassignedClusters,
ReassignClustersAlgorithm algorithm,
- Map<Cluster, Track> newMapShowerComponentToTrack,
- Map<Track, Cluster> tweakedTracksMatchedToClusters,
Map<Track, Set<Cluster>> newMapTrackToShowerComponents)
{
Set<Cluster> showerComponents = newMapJetToShowerComponents.get(jet);
- List<Cluster> reassignedClusters = reassignClustersToJet(jet, showerComponents, unassignedClusters, allSharedClusters, m_jetTolerance, algorithm, newMapShowerComponentToTrack, tweakedTracksMatchedToClusters,newMapTrackToShowerComponents);
+ List<Cluster> reassignedClusters = reassignClustersToJet(jet, showerComponents, unassignedClusters, allSharedClusters, m_jetTolerance, algorithm, newMapTrackToShowerComponents);
if (reassignedClusters != null && reassignedClusters.size()>0) {
for (Cluster clus : reassignedClusters) {
showerComponents.add(clus);
@@ -1268,14 +1289,14 @@
}
- protected List<Cluster> reassignClustersToTrack (Track tr, Collection<Cluster> initialShower, Collection<Cluster> unassignedClusters, List<SharedClusterGroup> allSharedClusters, double tolerance, ReassignClustersAlgorithm reassignAlgorithm, Map<Cluster, Track> newMapShowerComponentToTrack, Map<Track, Cluster> tweakedTracksMatchedToClusters, Map<Track, Set<Cluster>> newMapTrackToShowerComponents)
+ protected List<Cluster> reassignClustersToTrack (Track tr, Collection<Cluster> initialShower, Collection<Cluster> unassignedClusters, List<SharedClusterGroup> allSharedClusters, double tolerance, ReassignClustersAlgorithm reassignAlgorithm, Map<Track, Set<Cluster>> newMapTrackToShowerComponents)
{
Set<Track> tmpJet = new HashSet<Track>();
tmpJet.add(tr);
- return reassignClustersToJet(tmpJet, initialShower, unassignedClusters, allSharedClusters, tolerance, reassignAlgorithm, newMapShowerComponentToTrack, tweakedTracksMatchedToClusters, newMapTrackToShowerComponents);
+ return reassignClustersToJet(tmpJet, initialShower, unassignedClusters, allSharedClusters, tolerance, reassignAlgorithm, newMapTrackToShowerComponents);
}
- protected List<Cluster> reassignClustersToJet (Set<Track> jet, Collection<Cluster> initialShower, Collection<Cluster> unassignedClusters, List<SharedClusterGroup> allSharedClusters, double tolerance, ReassignClustersAlgorithm reassignAlgorithm, Map<Cluster, Track> newMapShowerComponentToTrack, Map<Track, Cluster> tweakedTracksMatchedToClusters, Map<Track, Set<Cluster>> newMapTrackToShowerComponents)
+ protected List<Cluster> reassignClustersToJet (Set<Track> jet, Collection<Cluster> initialShower, Collection<Cluster> unassignedClusters, List<SharedClusterGroup> allSharedClusters, double tolerance, ReassignClustersAlgorithm reassignAlgorithm, Map<Track, Set<Cluster>> newMapTrackToShowerComponents)
{
// Truth info debug
List<MCParticle> mcList = m_event.get(MCParticle.class, m_mcList);
@@ -1354,6 +1375,7 @@
break;
} else {
// Add this cluster and be willing to add more
+ if (m_debug) { System.out.println("DEBUG: Add this cluster and be willing to add more "+energyBefore+" -> "+energyAfter); }
previousIterationOK = true;
output.add(value);
}
@@ -1365,6 +1387,7 @@
output.add(value);
} else {
// Need more energy -- add cluster and keep going
+ if (m_debug) { System.out.println("DEBUG: Need more energy "+energyBefore+" -> "+energyAfter); }
output.add(value);
}
} else {
@@ -2055,6 +2078,101 @@
}
}
+ void scanForMuons(List<Track> trackList, Map<Track, BasicCluster> mipmap, List<Cluster> mipsMuon, List<Track> muonTracks){
+ for(Track tr : trackList){
+ List<MCParticle> m_truthList = getMCParticle(tr);
+ Hep3Vector track = (new BasicHep3Vector(tr.getMomentum()));
+ double p = track.magnitude();
+ double pT = Math.sqrt(track.x()*track.x() + track.y()*track.y());
+ HelixExtrapolationResult result = m_findCluster.performExtrapolation(tr);
+ Cluster mipcluster = mipmap.get(tr);
+
+ int layer = -1 ;
+ String lastsubdetName = "Not found";
+ Hep3Vector hitpos = new BasicHep3Vector();
+ Hep3Vector tangent = new BasicHep3Vector();
+ Hep3Vector tangentUnit = new BasicHep3Vector();
+
+ if(mipcluster == null || result == null) {
+ if(m_muonDebug){
+ if(mipcluster == null) aida.cloud1D("muon/mipcluster null").fill(pT);
+ if(result == null) aida.cloud1D("muon/result null").fill(pT);
+ }
+ if(pT > 3){//use IP to get tangent direction
+ tangent = track;
+ tangentUnit = VecOp.unit(tangent);
+ }else { continue; } // don't want to use this track to match with cluster in Muon Det.
+ }else if(mipcluster.getCalorimeterHits().size() > 0){
+ CalorimeterHit hit = mipcluster.getCalorimeterHits().get(mipcluster.getCalorimeterHits().size()-1);
+ layer = getLayer(hit);
+ lastsubdetName = hit.getSubdetector().getName();
+ hitpos = new BasicHep3Vector(hit.getPosition());
+ tangent = result.getTangent(hitpos);
+ tangentUnit = VecOp.unit(tangent);
+ } else {
+ System.out.println("Error: Has both mipcluster and extrapolation but no tangent");
+ continue;
+ }
+
+ if( mipsMuon != null){
+ double bestmatch = -1;
+ for(Cluster mipMuon : mipsMuon){
+ if( mipMuon.getCalorimeterHits().size() > 1){
+ CalorimeterHit muonhit1 = mipMuon.getCalorimeterHits().get(mipMuon.getCalorimeterHits().size()-1);
+ CalorimeterHit muonhit0 = mipMuon.getCalorimeterHits().get(0);
+ Hep3Vector muonpos1 = new BasicHep3Vector(muonhit1.getPosition());
+ Hep3Vector muonpos0 = new BasicHep3Vector(muonhit0.getPosition());
+ Hep3Vector trackMuon = VecOp.sub(muonpos1, muonpos0);
+ Hep3Vector trackMuonUnit = VecOp.unit(trackMuon);
+ double cos = VecOp.dot(tangentUnit, trackMuonUnit);
+ if(m_muonDebug){
+ if (cos == -1) {
+ System.out.println("tan: " + tangentUnit.x() + " " + tangentUnit.y() + " " + tangentUnit.z());
+ System.out.println("muon: " + trackMuonUnit.x() + " " + trackMuonUnit.y() + " " + trackMuonUnit.z());
+ }
+ }
+ if( cos > bestmatch) { bestmatch = cos; }
+ }
+ }
+ if(m_muonDebug){
+ aida.cloud1D("muon/match cos distribution",10000).fill(bestmatch);
+ }
+
+ double r = Math.sqrt(hitpos.x()*hitpos.x() + hitpos.y()*hitpos.y());
+ double z = hitpos.z();
+ if(m_muonDebug) {
+ System.out.println("Muon Debug: r= " + r + " z= " + z + " Det= " + lastsubdetName + " layer= " + layer);
+ }
+ //boolean tight = ( lastsubdetName == "MuonEndcap" || lastsubdetName == "MuonBarrel" );
+ //boolean loose = ( (layer >= 29 || layer <=33) && ( lastsubdetName.contains("HADEndcap") || lastsubdetName.contains("HADBarrel") )) || tight;
+ boolean match = bestmatch > 0.8;
+
+ for(MCParticle part : m_truthList){
+ if(part.getPDGID() == 13 || part.getPDGID() == -13) {
+ if(m_muonDebug){ aida.cloud1D("muon/MC muon").fill(part.getMomentum().magnitude());}
+ }
+ }
+
+ if(match){
+ muonTracks.add(tr);
+ boolean goodmatch = true;
+ if(m_muonDebug) {aida.cloud1D("muon/reconstructed muon").fill(p);}
+ for(MCParticle part : m_truthList){
+ if(part.getPDGID() == 13 || part.getPDGID() == -13) {
+ goodmatch = true;
+ break;
+ }
+ goodmatch = false;
+ }
+ if(m_muonDebug){
+ if(goodmatch) aida.cloud1D("muon/good matching muon").fill(p);
+ else {aida.cloud1D("muon/not good matching muon").fill(p);}
+ }
+ }
+ }
+ }
+ }
+
void scanForElectrons(List<Track> trackList, List<Cluster> photons, List<Cluster> electronClusters, List<Track> electronTracks) {
Map<Track,Cluster> electronCandidateMap = new HashMap<Track,Cluster>();
Set<Cluster> electronCandidateClusters = new HashSet<Cluster>();
@@ -2721,7 +2839,7 @@
List<Track> trackList, List<Track> tracksSortedByMomentum, List<Track> unmatchedTracksThatDontReachCalorimeter,
Map<Track, Track> mapOrigTrackToTweakedTrack,
Map<Track,Cluster> tracksMatchedToClusters, Map<Cluster, List<Track>> clustersMatchedToTracks,
- List<Track> electronTracks, List<Cluster> modifiedPhotonClusters, List<Cluster> electronClusters, List<Cluster> photonLikePhotons, List<Cluster> chargedHadronLikePhotons,
+ List<Track> electronTracks, List<Track> muonTracks, List<Cluster> modifiedPhotonClusters, List<Cluster> electronClusters, List<Cluster> photonLikePhotons, List<Cluster> chargedHadronLikePhotons,
List<Cluster> seedHadronLikePhotonClusters, List<Cluster> nonSeedHadronLikePhotonClusters, List<Cluster> nonSeedPhotonLikePhotonClusters,
Map<Track, Set<Cluster>> newMapTrackToShowerComponents, Map<Cluster, Track> newMapShowerComponentToTrack,
List<Cluster> linkableClusters,
@@ -2779,7 +2897,14 @@
ParticleType type_positron = ParticlePropertyManager.getParticlePropertyProvider().get(pdg_positron);
BaseParticleID pid_electron = new BaseParticleID(type_electron);
BaseParticleID pid_positron = new BaseParticleID(type_positron);
- double mass_electron = 0.0005; // rough!
+ double mass_electron = type_electron.getMass();
+ int pdg_muminus = 13;
+ int pdg_muplus = -13;
+ ParticleType type_muminus = ParticlePropertyManager.getParticlePropertyProvider().get(pdg_muminus);
+ ParticleType type_muplus = ParticlePropertyManager.getParticlePropertyProvider().get(pdg_muplus);
+ BaseParticleID pid_muminus = new BaseParticleID(type_muminus);
+ BaseParticleID pid_muplus= new BaseParticleID(type_muplus);
+ double mass_muon = type_muplus.getMass();
double windowEoverPveto = 2.5;
Set<Track> trackToTreatAsChargedParticles = new HashSet<Track>();
@@ -2882,6 +3007,7 @@
// FIXME: Break up MultipleTrackTrack tracks.
for (Track tr : tracksSortedByMomentum) {
boolean isElectron = electronTracks.contains(tr); // includes positrons
+ boolean isMuon = muonTracks.contains(tr);
// WARNING: Clusters may be entirely wrong for jets!
List<Cluster> clustersOfTrack = new Vector<Cluster>();
Set<Cluster> showerComponents = newMapTrackToShowerComponents.get(tr);
@@ -2917,6 +3043,9 @@
if (isElectron) {
currentParticleMass = mass_electron;
}
+ else if (isMuon) {
+ currentParticleMass = mass_muon;
+ }
double trackEnergySq = trackMomentumMag*trackMomentumMag + currentParticleMass*currentParticleMass;
HepLorentzVector fourMomentum = new BasicHepLorentzVector(Math.sqrt(trackEnergySq), trackMomentum);
part.set4Vector(fourMomentum);
@@ -2930,7 +3059,15 @@
part.addParticleID(pid_electron);
part.setParticleIdUsed(pid_electron);
}
- } else {
+ } else if(isMuon) {
+ if (part.getCharge()>0) {
+ part.addParticleID(pid_muplus);
+ part.setParticleIdUsed(pid_muplus);
+ } else {
+ part.addParticleID(pid_muminus);
+ part.setParticleIdUsed(pid_muminus);
+ }
+ } else {
if (part.getCharge()>0) {
part.addParticleID(pid_piplus);
part.setParticleIdUsed(pid_piplus);
