lcsim/src/org/lcsim/contrib/uiowa/MuonFinder
diff -N MuonFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonFinder.java 28 Sep 2008 06:32:06 -0000 1.1
@@ -0,0 +1,198 @@
+package org.lcsim.contrib.uiowa.MuonFinder;
+
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.util.*;
+import org.lcsim.event.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.base.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.recon.cluster.directedtree.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+import org.lcsim.recon.cluster.structural.likelihood.*;
+import org.lcsim.recon.cluster.structural.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ * Driver to find muon cluster and corresponding track.
+ *
+ * search mip cluster in Muon detector and assume that that cluster is from muon
+ * then find track matched with this mip cluster by looking at direction.
+ * Since pion can reach muon detector as well specially in endcap,
+ * it is required to find good quality mip cluster in muon detector.
+ */
+
+public class MuonFinder extends Driver{
+
+ protected String _outName;
+ protected String _tracklist;
+ protected String _outMap;
+ protected String _inMu;
+ protected String _inCal;
+ protected double _bestmatch = 0.8;
+ protected boolean _debug = false;
+
+ EventHeader m_event;
+
+ protected HelixExtrapolator m_findCluster = new org.lcsim.recon.pfa.identifier.LocalHelixExtrapolator();
+ protected boolean useMuonhits = false;
+
+ public MuonFinder(String tracklist, String inCal, String inMu, String outMap, String outName){
+ _inCal = inCal;
+ _inMu = inMu;
+ _outMap = outMap;
+ _outName = outName;
+ _tracklist = tracklist;
+ }
+
+ public void process(EventHeader event)
+ {
+ 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");
+ List<CalorimeterHit> allHitsHcalBarrel = event.get(CalorimeterHit.class, "HcalBarrDigiHits");
+ List<CalorimeterHit> allHitsHcalEndcap = event.get(CalorimeterHit.class, "HcalEndcapDigiHits");
+ calhits.addAll(allHitsEcalBarrel);
+ calhits.addAll(allHitsEcalEndcap);
+ calhits.addAll(allHitsHcalBarrel);
+ calhits.addAll(allHitsHcalEndcap);
+ List<CalorimeterHit> allHitsMUBarrel = event.get(CalorimeterHit.class, "MuonBarrDigiHits");
+ List<CalorimeterHit> allHitsMUEndcap = event.get(CalorimeterHit.class, "MuonEndcapDigiHits");
+ mudethits.addAll(allHitsMUBarrel);
+ mudethits.addAll(allHitsMUEndcap);
+ }else{
+ calhits.addAll(calHitMap.values());
+ mudethits.addAll(muHitMap.values());
+ }
+
+ m_findCluster.process(event); // picks up geometry
+
+ List<Track> trackList = event.get(Track.class, _tracklist);
+ List<Cluster> muonClusters = new Vector<Cluster>();
+ List<Track> muonTracks = new Vector<Track>();
+
+ //create mip cluster separately in Calorimetry and Muon detector.
+ //mip in CAL has track matched and muon mip in Muon detector is standalone.
+ //In the future, the muon mip will be replace with standalone muon track.
+ MipTrackMap finder = new MipTrackMap(m_findCluster);
+ Map<Track,BasicCluster> mipmap = finder.createCalMIPMap(calhits, trackList);
+ List<Cluster> muonmips = finder.createMIPMuDet(mudethits);
+
+ Map<Track,Set<Cluster>> outputmap = new HashMap<Track,Set<Cluster>>();
+
+ for(Cluster mumip : muonmips){
+ if( mumip.getCalorimeterHits().size() > 1){
+ CalorimeterHit muonhit1 = mumip.getCalorimeterHits().get(mumip.getCalorimeterHits().size()-1);
+ CalorimeterHit muonhit0 = mumip.getCalorimeterHits().get(0);
+ Hep3Vector muonpos1 = new BasicHep3Vector(muonhit1.getPosition());
+ Hep3Vector muonpos0 = new BasicHep3Vector(muonhit0.getPosition());
+ 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 bestmatch = -1;
+ Track bestTrack = null;
+ Cluster bestmipc = null;
+ for(Track tr : trackList){
+ 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 calmip = mipmap.get(tr);
+ Hep3Vector hitpos = new BasicHep3Vector();
+ Hep3Vector tangent = new BasicHep3Vector();
+ Hep3Vector tangentUnit = new BasicHep3Vector();
+
+ if(calmip == null || result == null) {
+ 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(calmip.getCalorimeterHits().size() > 0){
+ CalorimeterHit hit = calmip.getCalorimeterHits().get(calmip.getCalorimeterHits().size()-1);
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ String subdetName = hit.getSubdetector().getName();
+ //Followings are required to take pion track out for muon candidate
+ boolean muon1 = subdetName.contains("HADBarrel") || subdetName.contains("HADEndcap");
+ boolean muon2 = (layer > 29 && layer <40) ;
+ if(!(muon1 && muon2)) continue;
+ hitpos = new BasicHep3Vector(hit.getPosition());
+ tangent = result.getTangent(hitpos);
+ tangentUnit = VecOp.unit(tangent);
+ } else {
+ if(_debug) System.out.println("Error: Has both mipcluster and extrapolation but no tangent");
+ continue;
+ }
+ //look for angle between muon track and tangent vector
+ double cos = VecOp.dot(tangentUnit, trackMuonUnit);
+ //find best matched track
+ if( cos > bestmatch) {
+ bestmatch = cos;
+ bestTrack = tr;
+ bestmipc = calmip;
+ }
+ }
+ if(bestmatch > _bestmatch){
+ Set<Cluster> c = outputmap.get(bestTrack);
+ if (c == null) {
+ c = new HashSet<Cluster>();
+ outputmap.put(bestTrack, c);
+ //add if it has hits
+ if( bestmipc != null && bestmipc.getCalorimeterHits().size() > 0 ) c.add(bestmipc);
+ }
+ c.add(mumip);
+ }else { if(_debug) System.out.println("Warning: This cluster in MuDet has no track matched!");}
+ }else { if(_debug) System.out.println ("Warning: This muon mip has 0 number of hits!");}
+ }
+
+ List<CalorimeterHit> muonclustershits= new Vector<CalorimeterHit>();
+
+ for(Track tr : outputmap.keySet()){
+ muonTracks.add(tr);
+ BasicCluster tmpClus = new BasicCluster();
+ for(Cluster clus : outputmap.get(tr)){
+ for(CalorimeterHit hit : clus.getCalorimeterHits()){
+ tmpClus.addHit(hit);
+ muonclustershits.add(hit);
+ }
+ }
+ muonClusters.add(tmpClus);
+ }
+
+ HitMap muonMap = new HitMap(muonclustershits);
+ event.put(_outMap, muonMap);
+ int flag = 1<<LCIOConstants.CLBIT_HITS;
+ event.put(_outName, muonClusters, Cluster.class, flag );
+ event.put("MuonList", muonTracks);
+ }
+
+}
