lcsim/src/org/lcsim/contrib/uiowa/MuonFinder
diff -u -r1.5 -r1.6
--- MuonFinder.java 8 Oct 2008 20:06:41 -0000 1.5
+++ MuonFinder.java 11 Oct 2008 01:01:49 -0000 1.6
@@ -18,6 +18,8 @@
import org.lcsim.recon.cluster.mipfinder.*;
import org.lcsim.recon.cluster.clumpfinder.*;
import org.lcsim.geometry.IDDecoder;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.Detector;
import org.lcsim.util.lcio.LCIOConstants;
import org.lcsim.util.aida.AIDA;
/**
@@ -36,7 +38,7 @@
protected String _outMap;
protected String _inMu;
protected String _inCal;
- protected double _bestmatch = 0.9;
+ protected double _bestmatch = 0.8;
protected boolean _debug = false;
protected boolean _useFirstTwoLayer = true;
EventHeader m_event;
@@ -56,12 +58,14 @@
{
super.process(event);
m_event = event;
+
HitMap calHitMap = (HitMap) event.get(_inCal);
HitMap muHitMap = (HitMap) event.get(_inMu);
Set<CalorimeterHit> calhits = new HashSet<CalorimeterHit>();
Set<CalorimeterHit> mudethits = new HashSet<CalorimeterHit>();
+
if(useMuonhits){
List<CalorimeterHit> allHitsEcalBarrel = event.get(CalorimeterHit.class, "EcalBarrDigiHits");
List<CalorimeterHit> allHitsEcalEndcap = event.get(CalorimeterHit.class, "EcalEndcapDigiHits");
@@ -80,6 +84,10 @@
mudethits.addAll(muHitMap.values());
}
+ Set<Long> allHitIDs = new HashSet<Long>();
+ for (CalorimeterHit allhit : calhits) allHitIDs.add(allhit.getCellID());
+
+
m_findCluster.process(event); // picks up geometry
List<Track> trackList = event.get(Track.class, _tracklist);
@@ -96,7 +104,7 @@
Map<Track,Set<Cluster>> outputmap = new HashMap<Track,Set<Cluster>>();
if(_debug){
- System.out.println("muonmip size= " + muonmips.size()+" (should be 1 for signle sample)");
+ System.out.println("Muon mip size= " + muonmips.size()+" (should be 1 for signle sample)");
aida.cloud1D("muon/muon mip size").fill(muonmips.size());
}
for(Cluster mumip : muonmips){
@@ -144,8 +152,8 @@
muonpos1 = new BasicHep3Vector(hit1.getPosition());
}
- Hep3Vector trackMuon = VecOp.sub(muonpos1, muonpos0);
- Hep3Vector trackMuonUnit = VecOp.unit(trackMuon);
+ Hep3Vector Muon = VecOp.sub(muonpos1, muonpos0);
+ Hep3Vector MuonUnit = VecOp.unit(Muon);
double bestmatch = -1;
Track bestTrack = null;
@@ -158,12 +166,13 @@
HelixExtrapolationResult result = m_findCluster.performExtrapolation(tr);
Cluster calmip = mipmap.get(tr);
+ Hep3Vector last = new BasicHep3Vector();
+ Hep3Vector lastUnit = new BasicHep3Vector();
Hep3Vector hitpos = new BasicHep3Vector();
- Hep3Vector tangent = new BasicHep3Vector();
- Hep3Vector tangentUnit = new BasicHep3Vector();
- Hep3Vector firstlast = new BasicHep3Vector();
- Hep3Vector firstlastUnit = new BasicHep3Vector();
-
+ Hep3Vector lastpoint = new BasicHep3Vector();
+
+ boolean muon = false;
+ int isolated = 0;
if(calmip == null || result == null) {
if(_debug) System.out.println("Null mip or extrapolation");
continue;
@@ -174,30 +183,92 @@
int layer = id.getLayer();
String subdetName = hit.getSubdetector().getName();
//Followings are required to take pion track out for muon candidate
- boolean muon = subdetName.contains("HADBarrel") || (subdetName.contains("HADEndcap") && (layer > 29 && layer <40) );
+ muon = subdetName.contains("HADBarrel") || (subdetName.contains("HADEndcap") && (layer > 29) );
hitpos = new BasicHep3Vector(hit.getPosition());
- if(!(muon)) continue;
- tangent = result.getTangent(hitpos);
- tangentUnit = VecOp.unit(tangent);
- firstlast = VecOp.sub(muonpos0, hitpos);
- firstlastUnit = VecOp.unit(firstlast);
+ //Using extrapolated track
+ int lastLayerBarrel = getLastLayer("HADBarrel");
+ int lastLayerEndcap = getLastLayer("HADEndcap");
+ if(lastLayerBarrel != lastLayerEndcap) {System.out.println("ERROR : number of layer different");}
+ int endingLayer = lastLayerBarrel;
+ boolean isBarrel = true;
+ for(int i=0 ; i < 20 ;i++){
+ endingLayer = endingLayer - 1;
+ IDDecoder tid = null;
+ Hep3Vector tpoint = new BasicHep3Vector();
+ Hep3Vector tBarrel = null;
+ if(isBarrel) {
+ tBarrel = result.extendToHCALBarrelLayer(endingLayer);
+ }
+ Hep3Vector tEndcap = result.extendToHCALEndcapLayer(endingLayer);
+ if(tBarrel != null && tEndcap == null){
+ tpoint = tBarrel;
+ }else if(tBarrel == null && tEndcap != null){
+ tpoint = tEndcap;
+ isBarrel = false;
+ }else if(tBarrel != null && tEndcap != null){
+ //System.out.println("Both have it");
+ tpoint = tBarrel;
+ }else{
+ //System.out.println("Both don't have it");
+ for(int l=0; l<endingLayer+1 ; l++){
+ int nextlayer = endingLayer-l;
+ if(isBarrel) {tpoint = result.extendToHCALBarrelLayer(nextlayer);}
+ else {tpoint = result.extendToHCALEndcapLayer(nextlayer);}
+ if(tpoint != null){
+ endingLayer = nextlayer;
+ break;
+ }
+ }
+ }
+
+ if(isBarrel){tid = event.getDetector().getDecoder("HcalBarrHits");}
+ else {tid = event.getDetector().getDecoder("HcalEndcapHits");}
+ long cellID = tid.findCellContainingXYZ(tpoint);
+ tid.setID(cellID);
+ long[] neighbours = tid.getNeighbourIDs(0,5,5);
+ int count = 0;
+ if(calHitMap.containsKey(cellID)) {count++;}
+ for(long neighbour : neighbours){
+ if(calHitMap.containsKey(neighbour)){count++;}
+ }
+
+ if(i==0){
+ lastpoint = tpoint;
+ last = result.getTangent(tpoint);
+ lastUnit = VecOp.unit(last);
+ }
+ double rpos = Math.sqrt(tpoint.x()*tpoint.x() + tpoint.y()*tpoint.y());
+ System.out.println("Extrapolated track at " + endingLayer + " : r= "+rpos+" z= "+tpoint.z() + " hit= " + count );
+ if(endingLayer == 0){
+ endingLayer = 40;
+ isBarrel = false;
+ }
+ if(count > 0 && count < 3) isolated++;
+ }
+
} else {
if(_debug) System.out.println("Error: Has both mipcluster and extrapolation but no tangent");
continue;
}
- double cos0 = VecOp.dot(tangentUnit, trackMuonUnit);
- double cos1 = VecOp.dot(tangentUnit, firstlastUnit);
- double cos2 = VecOp.dot(trackMuonUnit, firstlastUnit);
+ //require last position or 9 isolated layer in last 20 layers obtained by extrapolated track.
+ if(!(muon || (isolated > 9) )) continue;
+ Hep3Vector link = VecOp.sub(muonpos0, lastpoint);
+ Hep3Vector linkUnit = VecOp.unit(link);
+ double cos0 = VecOp.dot(lastUnit, MuonUnit);
+ double cos1 = VecOp.dot(lastUnit, linkUnit);
+ double cos2 = VecOp.dot(MuonUnit, linkUnit);
double cos = (cos0 + cos1 + cos2)/3;
//find best matched track
+ //compare the combination of three directions
if( cos > bestmatch ) {
if(_debug){
- aida.cloud1D("muon/cos_tangent_muon").fill(cos0);
- aida.cloud1D("muon/cos_tangent_firstlast").fill(cos1);
- aida.cloud1D("muon/cos_muon_firstlast").fill(cos2);
+ aida.cloud1D("muon/cos_extrapolatedtrack").fill(cos0);
+ aida.cloud1D("muon/cos_extrapolatedtrack_link").fill(cos1);
+ aida.cloud1D("muon/cos_muon_link").fill(cos2);
aida.cloud1D("muon/combined cos").fill(cos);
System.out.println("cos0= " + cos0 + " cos1= " + cos1 + " cos2= " + cos2);
+ System.out.println("average cos= " + cos);
}
bestmatch = cos;
bestTrack = tr;
@@ -257,5 +328,11 @@
return meanPosition;
}
+
+ public int getLastLayer(String s){
+ Detector det = m_event.getDetector();
+ CylindricalCalorimeter detname = ((CylindricalCalorimeter) det.getSubdetectors().get(s));
+ return detname.getLayering().getLayerCount();
+ }
}
