lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.50 -r1.51
--- ReclusterDTreeDriver.java 22 Sep 2008 23:59:25 -0000 1.50
+++ ReclusterDTreeDriver.java 23 Sep 2008 16:45:40 -0000 1.51
@@ -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.50 2008/09/22 23:59:25 tjkim Exp $
+ * @version $Id: ReclusterDTreeDriver.java,v 1.51 2008/09/23 16:45:40 tjkim Exp $
* @author Mat Charles <[log in to unmask]>
*/
@@ -284,6 +284,7 @@
List<Cluster> mipsMuon = new Vector<Cluster>();
mipsMuon.addAll(event.get(Cluster.class, "OldMipsInsideTreesMCALforMuonID"));
mipsMuon.addAll(event.get(Cluster.class, "NewMipsInsideTreesMCALforMuonID"));
+ mipsMuon.addAll(event.get(Cluster.class, "ClumpsInsideTreesMCALforMuonID"));
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.");
}
@@ -2080,10 +2081,12 @@
void scanForMuons(List<Track> trackList, Map<Track, BasicCluster> mipmap, List<Cluster> mipsMuon, List<Track> muonTracks){
for(Track tr : trackList){
+ if(m_muonDebug){ System.out.println("start scanForMuons for this track");}
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());
+ double costheta = Math.abs(track.z()/p);
HelixExtrapolationResult result = m_findCluster.performExtrapolation(tr);
Cluster mipcluster = mipmap.get(tr);
@@ -2097,10 +2100,17 @@
if(m_muonDebug){
if(mipcluster == null) aida.cloud1D("muon/mipcluster null").fill(pT);
if(result == null) aida.cloud1D("muon/result null").fill(pT);
+ if(result == null && mipcluster == null) aida.cloud1D("muon/both bull").fill(pT);
}
- if(pT > 3){//use IP to get tangent direction
- tangent = track;
- tangentUnit = VecOp.unit(tangent);
+ if(pT > 1){//use intercept point or IP to get tangent direction
+ if( result != null){
+ Hep3Vector interceptPoint = result.getInterceptPoint();
+ tangent = result.getTangent(interceptPoint);
+ tangentUnit = VecOp.unit(tangent);
+ }else{
+ 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);
@@ -2122,17 +2132,23 @@
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 trackMuon = new BasicHep3Vector();
+ if(muonpos0.magnitude() > muonpos1.magnitude()){
+ trackMuon = VecOp.sub(muonpos0, muonpos1);
+ }else{
+ 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 < 0) {//how come direction is opposite?
+ System.out.println("tan: "+tangentUnit.x()+" "+tangentUnit.y()+" "+tangentUnit.z());
+ System.out.println("muon: "+trackMuonUnit.x()+" "+trackMuonUnit.y()+" "+trackMuonUnit.z());
+ System.out.println("cos: " + cos );
}
}
if( cos > bestmatch) { bestmatch = cos; }
- }
+ }else { System.out.println ("This muon mip has 0 number of hits!");}
}
if(m_muonDebug){
aida.cloud1D("muon/match cos distribution",10000).fill(bestmatch);
@@ -2142,21 +2158,29 @@
double z = hitpos.z();
if(m_muonDebug) {
System.out.println("Muon Debug: r= " + r + " z= " + z + " Det= " + lastsubdetName + " layer= " + layer);
+ System.out.println("Cos theta= " + bestmatch);
}
- //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(m_muonDebug){
+ aida.cloud1D("muon/FastTrack muon").fill(p);
+ aida.cloud1D("muon/FastTrack pT muon").fill(pT);
+ aida.cloud1D("muon/FastTrack costheta").fill(costheta);
+ }
}
}
if(match){
muonTracks.add(tr);
boolean goodmatch = true;
- if(m_muonDebug) {aida.cloud1D("muon/reconstructed muon").fill(p);}
+ if(m_muonDebug) {
+ aida.cloud1D("muon/reconstructed muon").fill(p);
+ aida.cloud1D("muon/reconstructed pT muon").fill(pT);
+ aida.cloud1D("muon/reconstructed costheta").fill(costheta);
+ }
for(MCParticle part : m_truthList){
if(part.getPDGID() == 13 || part.getPDGID() == -13) {
goodmatch = true;
@@ -2165,11 +2189,19 @@
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);}
+ if(goodmatch) {
+ aida.cloud1D("muon/good matching muon").fill(p);
+ aida.cloud1D("muon/good matching pT muon").fill(pT);
+ aida.cloud1D("muon/good matching costheta").fill(costheta);
+ }
+ else {
+ aida.cloud1D("muon/not good matching muon").fill(p);
+ aida.cloud1D("muon/not good matching pT muon").fill(pT);
+ aida.cloud1D("muon/not good matching costheta").fill(costheta);
+ }
}
}
- }
+ }else { System.out.println ("No muon mip cluster found : Null");}
}
}