lcsim/src/org/lcsim/contrib/uiowa
diff -N MuonStandAlone.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonStandAlone.java 24 Nov 2008 22:29:14 -0000 1.1
@@ -0,0 +1,38 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.loop.LCIODriver;
+import org.lcsim.util.*;
+import org.lcsim.event.*;
+import org.lcsim.recon.cluster.structural.likelihood.LikelihoodEvaluatorWrapper;
+import org.lcsim.recon.cluster.directedtree.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.recon.pfa.output.EnergySumPlotter;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.digisim.*;
+import org.lcsim.mc.fast.tracking.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.recon.pfa.identifier.*;
+
+public class MuonStandAlone extends Driver
+{
+ public MuonStandAlone() {
+ // Prepare to run PFA: Tracks (includes DigiSim)
+ add(new org.lcsim.digisim.DigiPackageDriver());
+ add(new org.lcsim.recon.tracking.seedtracker.ReconTracking.SiD02ReconTrackingDriver());
+
+ // Run MuonFinderDriver
+ org.lcsim.contrib.uiowa.MuonFinder.MuonFinderDriver muonFinderDriver = new org.lcsim.contrib.uiowa.MuonFinder.MuonFinderDriver("Tracks", "Muons");
+ add(muonFinderDriver);
+ add(new org.lcsim.util.loop.LCIODriver("rerun.slcio"));
+ }
+
+ int count = 0;
+ public void process(EventHeader event) {
+ System.out.println("DEBUG: Looking at event "+count); count++;
+ super.process(event);
+ }
+}
lcsim/src/org/lcsim/contrib/uiowa/MuonFinder
diff -N MuonFinderDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonFinderDriver.java 24 Nov 2008 22:29:14 -0000 1.1
@@ -0,0 +1,349 @@
+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.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.util.aida.AIDA;
+/**
+ * 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 MuonFinderDriver extends Driver{
+ protected AIDA aida = AIDA.defaultInstance();
+ protected String _tracklist;
+ protected String _muonlist;
+ protected double _bestmatch = 0.8;
+ protected boolean _debug = false;
+ protected boolean _useFirstTwoLayer = true;
+ EventHeader m_event;
+
+ protected HelixExtrapolator _findCluster = new org.lcsim.recon.pfa.identifier.TrackHelixPlusHitExtrapolator();
+
+ public MuonFinderDriver(String tracklist, String muonlist){
+ _tracklist = tracklist;
+ _muonlist = muonlist;
+ }
+
+ public void process(EventHeader event)
+ {
+ super.process(event);
+ m_event = event;
+
+ Set<CalorimeterHit> calhits = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> mudethits = new HashSet<CalorimeterHit>();
+
+ 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, "CorrMuonBarrDigiHits");
+ List<CalorimeterHit> allHitsMUEndcap = event.get(CalorimeterHit.class, "CorrMuonEndcapDigiHits");
+ mudethits.addAll(allHitsMUBarrel);
+ mudethits.addAll(allHitsMUEndcap);
+
+ HitMap calHitMap = new HitMap();
+ HitMap muHitMap = new HitMap();
+ for (CalorimeterHit hit : calhits) {
+ long id = hit.getCellID();
+ calHitMap.put(id, hit);
+ }
+
+ for (CalorimeterHit hit : mudethits) {
+ long id = hit.getCellID();
+ muHitMap.put(id, hit);
+ }
+
+ Set<Long> allHitIDs = new HashSet<Long>();
+ for (CalorimeterHit allhit : calhits) allHitIDs.add(allhit.getCellID());
+
+ _findCluster.process(event); // picks up geometry
+
+ List<Track> trackList = event.get(Track.class, _tracklist);
+
+ //Create mip cluster separately in Calorimetry and Muon detector.
+ //Mip in CAL has track matched
+ //But Mip in Muon detector is standalone.
+ //In the future, the muon Mip will be replace with standalone muon track.
+ MipTrackMap finder = new MipTrackMap(_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>>();
+ if(_debug){
+ 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){
+ if( mumip.getCalorimeterHits().size() > 1){
+ Hep3Vector muonpos0 = new BasicHep3Vector();
+ Hep3Vector muonpos1 = new BasicHep3Vector();
+ //Muon detector has 20cm iron between 2 layers. The direction is not affected by the segmentation.
+ //Trying to use first two layers to identify direction intead of using first and last hit.
+ //In endcap, track curves so it will bring wrong direction if we use first and last hit.
+ //First two layers will give better direction.
+ if(_useFirstTwoLayer){
+ SortedMap<Integer, Set<CalorimeterHit>> exam = new TreeMap();
+ int prelayer = -1 ;
+ int i = -1;
+ String predetName = null;
+ for(CalorimeterHit h : mumip.getCalorimeterHits()){
+ IDDecoder id = h.getIDDecoder();
+ id.setID(h.getCellID());
+ int layer = id.getLayer();
+ if(prelayer != layer || prelayer == -1){
+ predetName = h.getSubdetector().getName();
+ i++;
+ if(i==2) break;
+ }
+
+ Set<CalorimeterHit> hits = exam.get(layer);
+ if(hits == null){
+ hits = new HashSet<CalorimeterHit>();
+ exam.put(layer,hits);
+ }
+ String subdetName = h.getSubdetector().getName();
+ if(predetName.contains(subdetName)) hits.add(h);
+ prelayer = layer;
+ }
+
+ Collection<CalorimeterHit> hits0 = exam.get(exam.firstKey());
+ Collection<CalorimeterHit> hits1 = exam.get(exam.lastKey());
+ muonpos0 = getHitMeanPosition(hits0);
+ muonpos1 = getHitMeanPosition(hits1);
+
+ }else{
+ CalorimeterHit hit0 = mumip.getCalorimeterHits().get(0);
+ CalorimeterHit hit1 = mumip.getCalorimeterHits().get(mumip.getCalorimeterHits().size()-1);
+ muonpos0 = new BasicHep3Vector(hit0.getPosition());
+ muonpos1 = new BasicHep3Vector(hit1.getPosition());
+ }
+
+ Hep3Vector Muon = VecOp.sub(muonpos1, muonpos0);
+ Hep3Vector MuonUnit = VecOp.unit(Muon);
+
+ 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 = _findCluster.performExtrapolation(tr);
+ Cluster calmip = mipmap.get(tr);
+ Hep3Vector last = new BasicHep3Vector();
+ Hep3Vector lastUnit = new BasicHep3Vector();
+ Hep3Vector hitpos = 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;
+ }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
+ muon = subdetName.contains("HADBarrel") || (subdetName.contains("HADEndcap") && (layer > 29) );
+ hitpos = new BasicHep3Vector(hit.getPosition());
+ //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;
+ boolean hasLastPoint = 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(tpoint == null) {
+ if(_debug) System.out.println("Null extrapolated track point");
+ hasLastPoint = false;
+ 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());
+ if(_debug) 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++;
+ }
+ if(!hasLastPoint) continue;
+ } else {
+ if(_debug) System.out.println("Error: Has both mipcluster and extrapolation but no tangent");
+ continue;
+ }
+
+ //Require last position or 7 isolated layer in last 10 layers obtained by extrapolated track.
+ //This helps remove pion which is showering in CAL and also take the MIP which is contaminated by
+ //other shower so it stopped in the middle of CAL.
+ if(!(muon || (isolated > 6) )) continue;
+ if (!muon) {
+ // Skip these cases to avoid pion contamination.
+ if(_debug) System.out.println("Mip doesn't reach in last 10 layers");
+ 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_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;
+ 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!");
+ System.out.println("best match= " + bestmatch);
+ }
+ }
+ }else { if(_debug) System.out.println ("Warning: This muon mip has 0 number of hits!");}
+ }
+
+ List<CalorimeterHit> muonclustershits= new Vector<CalorimeterHit>();
+ Map<Track, Cluster> MuonMapTrackToCluster = new HashMap<Track,Cluster>();
+
+ for(Track tr : outputmap.keySet()){
+ BasicCluster tmpClus = new BasicCluster();
+ BasicCluster muonClus = new BasicCluster();
+ Cluster calmip = mipmap.get(tr);
+ for(Cluster clus : outputmap.get(tr)){
+ for(CalorimeterHit hit : clus.getCalorimeterHits()){
+ if(calmip != clus) { muonClus.addHit(hit); }
+ tmpClus.addHit(hit);
+ muonclustershits.add(hit);
+ }
+ }
+
+ //Check Quality for mip in Muon Detector only
+ SimpleMipQualityDecision check = new SimpleMipQualityDecision();
+ boolean passesCheck = check.valid(muonClus);
+ if(passesCheck) {
+ MuonMapTrackToCluster.put(tr,tmpClus);
+ }
+ }
+
+ HitMap muonMap = new HitMap(muonclustershits);
+
+ List<Track> muonTracks = new Vector<Track>(MuonMapTrackToCluster.keySet());
+ int flag = 1<<LCIOConstants.CLBIT_HITS;
+ event.put(_muonlist, muonTracks, Track.class, flag);
+ }
+
+ public Hep3Vector getHitMeanPosition(Collection<CalorimeterHit> hits){
+ Hep3Vector hitPosition = new BasicHep3Vector();
+ double size = hits.size();
+ double norm = 1/size;
+ for(CalorimeterHit hit : hits) {
+ hitPosition = VecOp.add(hitPosition, new BasicHep3Vector(hit.getPosition()));
+ }
+ Hep3Vector meanPosition = VecOp.mult(norm, hitPosition);
+
+ return meanPosition;
+ }
+
+ public int getLastLayer(String s){
+ Detector det = m_event.getDetector();
+ CylindricalCalorimeter detname = ((CylindricalCalorimeter) det.getSubdetectors().get(s));
+ return detname.getLayering().getLayerCount();
+ }
+
+}