Commit in lcsim/src/org/lcsim/recon/cluster/muonfinder on MAIN
MipTrackMap.java+523added 1.1
MuonFinder.java+358added 1.1
MuonFinderWrapper.java+109added 1.1
SimpleMipQualityDecision.java+78added 1.1
+1068
4 added files
MJC: Propagate MuonFinder to stable

lcsim/src/org/lcsim/recon/cluster/muonfinder
MipTrackMap.java added at 1.1
diff -N MipTrackMap.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MipTrackMap.java	22 Oct 2008 23:24:34 -0000	1.1
@@ -0,0 +1,523 @@
+package org.lcsim.recon.cluster.muonfinder;
+
+import java.util.*; 
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import hep.physics.vec.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.base.*;
+import org.lcsim.recon.cluster.structural.likelihood.*;
+import org.lcsim.recon.cluster.structural.*;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * MIP finder.
+ * make map for track and mip in CAL
+ * find sub mips in Muon Detector with quality cut.
+ * later need to make it as one cluster if it is assocated with the same track. 
+ *
+ * @author [log in to unmask]
+ * @version $Id: MipTrackMap.java,v 1.1 2008/10/22 23:24:34 mcharles Exp $
+ */
+
+
+public class MipTrackMap extends Driver {
+    protected Set<CalorimeterHit> m_allhits;
+    protected HelixExtrapolator m_extrap;
+    protected int exam = 0; //to test neighbour
+    protected boolean debug = false;
+    protected TrackClusterMatcher m_trackClusterMatcher;
+    protected HitMap hitsInTree = null;
+
+    protected double m_newMipFinderRadius = 100.0;
+    protected int m_minHitsToBeTreatedAsCluster = 5;
+
+    protected boolean m_removePoorQualityMips = true;
+    protected boolean m_doOldMipsOutsideTrees = false;
+    protected int m_dU =2;
+    protected int m_dV =2;
+    protected int m_dL =1;
+
+    protected boolean use_DTree = false;
+
+    public MipTrackMap(HelixExtrapolator extrap) {
+		m_extrap = extrap;
+                initTrackMatch();
+    }
+	 
+    public Map<Track, BasicCluster> createCalMIPMap(Set<CalorimeterHit> allhits, List<Track> tracklist) {
+        m_allhits = allhits;
+        List<Track> trackList = tracklist;
+        List<Cluster> allMatchableClusters = createClusters(m_allhits);
+	Map<Track, BasicCluster> MapTrkToMIP = new HashMap<Track, BasicCluster>();
+	for(Track tr : trackList){
+            exam = 0;
+            Hep3Vector track = (new BasicHep3Vector(tr.getMomentum()));
+            double mo = track.magnitude();
+            Cluster matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, allMatchableClusters);
+            BasicCluster mipclus = new BasicCluster();  //for new mip cluster   
+	    HelixExtrapolationResult result = m_extrap.performExtrapolation(tr);
+	    Hep3Vector interceptPoint = null; // for debugging
+	    if (result != null) {
+		interceptPoint = result.getInterceptPoint();
+	    }
+	    Cluster seed = matchedCluster;
+
+	    //It might be needed to use intecept point instead of using first hit of seed
+	    //so that we can find the exact first hit corresponding to track
+        
+            if(seed == null) continue; 
+	    // first hit of seed 
+	    mipclus.addHit(seed.getCalorimeterHits().get(0));
+	    //test for neighbour
+            examPosition(seed.getCalorimeterHits().get(0));
+
+	    int lastIterationWithFoundHit = -1;
+	    int maxSkippedLayers = 3;
+	    for( int iIteration = 0 ; iIteration < 71 ; iIteration++){ // for sid02 calorimeter. 31 for Ecal 40 for Hacl.
+		//To get the information from last hit
+		CalorimeterHit hit = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-1);
+		IDDecoder id = hit.getIDDecoder();
+		id.setID(hit.getCellID());
+		double x = id.getPosition()[0];
+		double y = id.getPosition()[1];
+		double z = id.getPosition()[2];
+		double p = Math.sqrt(x*x+y*y+z*z);			
+		double r = Math.sqrt(x*x+y*y);
+		double lastlayer = id.getLayer();
+		String lastsubdetName = hit.getSubdetector().getName();
+		// Check for neighbouring hits in a 5x5x(2n+1) grid
+		// where n is the number of iterations since we last saw a hit.
+		// For example, if we saw a hit in the previous layer/iteration
+		// it'll be a 5x5x3 grid.
+		Set<Long> nearIDs = new HashSet<Long>();
+		int countIterationsSinceLastHit = (iIteration - lastIterationWithFoundHit);
+		if (countIterationsSinceLastHit <= 0) { throw new AssertionError("Book-keeping error"); }
+		if (countIterationsSinceLastHit > maxSkippedLayers+1) { throw new AssertionError("Book-keeping error"); }
+		long[] nearbyHitArray = id.getNeighbourIDs(countIterationsSinceLastHit, 3, 3);
+		for (long lastnearID : nearbyHitArray ) { nearIDs.add(lastnearID); }
+
+		Hep3Vector last0 = new BasicHep3Vector(hit.getPosition());
+		Hep3Vector last1 = new BasicHep3Vector();
+		if( mipclus.getCalorimeterHits().size() > 1){
+		    CalorimeterHit hit1 = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-2);
+	            last1 = new BasicHep3Vector(hit1.getPosition());
+		}
+
+		//Searh cur hit in cur layer and arrange in angular order
+		SortedMap<Double, CalorimeterHit> sortedHitbyPos= new TreeMap();
+		for (CalorimeterHit curhit : m_allhits){
+		    IDDecoder curid = curhit.getIDDecoder();
+		    curid.setID(curhit.getCellID());
+		    double curx = curid.getPosition()[0];
+		    double cury = curid.getPosition()[1];
+		    double curz = curid.getPosition()[2];
+		    double curp = Math.sqrt(curx*curx+cury*cury+curz*curz);
+		    double curr = Math.sqrt(curx*curx+cury*cury);  
+		    Hep3Vector curpos = new BasicHep3Vector(curx,cury,curz);
+		    int curlayer = curid.getLayer();
+		    String cursubdetName = curhit.getSubdetector().getName();
+		    //use geometric information and layer number to find cur hit
+		    boolean InSide = (curlayer > lastlayer) && cursubdetName.contains(lastsubdetName) && nearIDs.contains(curhit.getCellID());
+                    boolean EcalToHcalBarrel= (curid.getLayer()==0 && cursubdetName.contains("HADBarrel")) && lastsubdetName.contains("EMBarrel"); // ECal Barrel -> HCal Barrel
+                    boolean EcalToHcalEndcap= (curid.getLayer()==0 && cursubdetName.contains("HADEndcap")) && lastsubdetName.contains("EMEndcap"); // ECal Endcap -> HCal Endcap
+                    boolean BarrelEcalToEndHcal = (curid.getLayer()==0 && cursubdetName.contains("HADEndcap")) && lastsubdetName.contains("EMBarrel"); // ECal Barrel -> HCal Endcap (This is no way to HCal Barrel directly from ECal Endcap
+		    boolean EndToBarrelEcal = (curid.getLayer()==0 && cursubdetName.contains("EMBarrel")) && lastsubdetName.contains("EMEndcap"); // ECal Endcap -> ECal Barrel 
+                    boolean EndToBarrelHcal = (curid.getLayer()==0 && cursubdetName.contains("HADBarrel")) && lastsubdetName.contains("HADEndcap"); // HCal Endcap -> HCal Barrel
+		    if(EcalToHcalBarrel || EcalToHcalEndcap || InSide || EndToBarrelEcal || EndToBarrelHcal || BarrelEcalToEndHcal){
+                        Hep3Vector curFromTangent = result.getTangent(curpos);
+			Hep3Vector curFromLastTwo = VecOp.sub(curpos, last0);
+                        Hep3Vector cur = null;
+                        if(curFromTangent != null) {
+                            cur = curFromTangent;
+                        }else {
+                            cur = curFromLastTwo;  
+                            if (debug) { System.out.println("No getTangent From current"); }
+                        }
+			Hep3Vector curUnit = VecOp.unit(cur);
+
+                        Hep3Vector lastFromTangent = result.getTangent(last0);
+			Hep3Vector lastFromLastTwo = null;
+			if (mipclus.getCalorimeterHits().size() > 1) { 
+			    lastFromLastTwo = VecOp.sub(last0, last1);
+			} else {
+			    //Since there is only one hit we have at the first, we use a extrapolrated tangent vector 	
+			    if (result != null) {
+				lastFromLastTwo = result.getTangent();
+			    }
+			    if (lastFromLastTwo == null) {
+			//	 No extrapolation info and only one hit -- fall back to using
+			//	 a straight line from the origin.
+				lastFromLastTwo = VecOp.sub(last0, new BasicHep3Vector(0,0,0));
+                                 
+			    }
+			}
+                        Hep3Vector last = null;                       
+                        if(lastFromTangent != null) {
+                            last = lastFromTangent;
+                        }else {
+                            last = lastFromLastTwo;
+                            if (debug) { System.out.println("No getTangent From last"); }
+                        }
+
+			Hep3Vector lastUnit = VecOp.unit(last);
+
+                        double dot = VecOp.dot(lastUnit, curUnit);
+                        if(dot > 1) dot = 1.0;
+			double a = Math.acos(dot);
+			double d = VecOp.sub(curpos,last0).magnitude();
+			if(a < 1.05 && d < 100){ //next hit in the range of angle 90 degree and less than 100mm(?)
+			    sortedHitbyPos.put(a, curhit);
+			}
+		    }		
+		}		
+		
+		if(sortedHitbyPos.size() == 0) {
+		    // No hit found in this layer.
+		    if (countIterationsSinceLastHit <= maxSkippedLayers) {
+			continue; // go to next layer 
+		    } else {
+			break; // stop -- ran out of hits and haven't seen any for several iterations
+		    }
+		} else {
+		    // Found one or more hits
+		    lastIterationWithFoundHit = iIteration;
+		}
+		CalorimeterHit hitAdd = sortedHitbyPos.get(sortedHitbyPos.firstKey()); //add the hit with a smallest angle 
+		examPosition(hitAdd);
+		mipclus.addHit(hitAdd);
+                
+		if (exam > 4) {
+		    if(debug) System.out.println("Debug: showering is already started " + exam + " layers ago" );
+                    double numRemove = 0;
+		    if(mipclus.getCalorimeterHits().size() < 6) {
+                        numRemove = mipclus.getCalorimeterHits().size() - 1 ; 
+                    } else numRemove = 5;
+		    for(int k=0; k < numRemove ; k++) {
+			CalorimeterHit toRemove = mipclus.getCalorimeterHits().get(mipclus.getCalorimeterHits().size()-1);
+			mipclus.removeHit(toRemove);
+		    }
+		    break;
+		}
+	    } //layer loop after mip
+
+	    if(debug){
+		int k=0;
+		if(mipclus.getCalorimeterHits().size() < 2) k = mipclus.getCalorimeterHits().size();
+		else k = 2;
+		for(int i =0; i < k ;i++){
+		    CalorimeterHit hit = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-1-i);
+		    Hep3Vector v = new BasicHep3Vector(hit.getPosition());
+		    System.out.println("Debug: last "+i+" pos= "+v.magnitude()+" (" +v.x() + " " + v.y() + " " + v.z()+ ")" );
+		}
+		System.out.println("cal MIP size= " + mipclus.getCalorimeterHits().size());
+	    }		
+	    MapTrkToMIP.put(tr, mipclus);
+	} //track loop
+
+	return MapTrkToMIP;
+    }
+
+    public List<Cluster> createMIPMuDet(Set<CalorimeterHit> muhits){
+        removeUpperSensitiveLayer(muhits);
+        hitsInTree = new HitMap(muhits);
+        List<Cluster> muonmip = new Vector<Cluster>();
+
+        if(!use_DTree){
+            Set<Long> allHitIDs = new HashSet<Long>();
+            Set<CalorimeterHit> allhits = new HashSet<CalorimeterHit>(); 
+            Set<CalorimeterHit> removedhits = new HashSet<CalorimeterHit>();
+            for(CalorimeterHit hit : muhits){
+                if(removedhits.contains(hit)) continue;
+                SortedMap<Double, Set<CalorimeterHit>> sortedHitbyPos= new TreeMap();
+                IDDecoder id = hit.getIDDecoder();
+                double layer = id.getLayer();
+                Hep3Vector hitPos = new BasicHep3Vector(hit.getPosition());
+                Hep3Vector hitPosUnit = VecOp.unit(hitPos);
+
+                Set<CalorimeterHit> hits = sortedHitbyPos.get(hitPos.magnitude()); 
+                if(hits == null){
+                    hits = new HashSet<CalorimeterHit>();
+                    sortedHitbyPos.put(hitPos.magnitude(),hits);
+                }
+                hits.add(hit);
+
+                removedhits.add(hit);
+                SortedMap<Double, CalorimeterHit> sortedHitbyAng= new TreeMap();
+                for(CalorimeterHit subhit : muhits){
+                    if(removedhits.contains(subhit)) continue;
+                    IDDecoder subid = subhit.getIDDecoder();
+                    double sublayer = subid.getLayer();
+                    Hep3Vector subhitPos = new BasicHep3Vector(subhit.getPosition()); 
+                    Hep3Vector subhitPosUnit = VecOp.unit(subhitPos);
+                    double cos = VecOp.dot(hitPosUnit, subhitPosUnit);
+                    sortedHitbyAng.put(cos,subhit);
+                    if(cos > 0.96) {
+
+                        Set<CalorimeterHit> subhits = sortedHitbyPos.get(subhitPos.magnitude());
+                        if(subhits == null){
+                            subhits = new HashSet<CalorimeterHit>();
+                            sortedHitbyPos.put(subhitPos.magnitude(),subhits);
+                        }
+                        if(subhit == null) { System.out.println("ERROR no subhit"); }
+                        subhits.add(subhit);
+                        removedhits.add(subhit);
+                    }
+                } 
+
+                BasicCluster tmpClus = new BasicCluster();
+                for(Set<CalorimeterHit> c : sortedHitbyPos.values()){
+                    for(CalorimeterHit h : c){
+                        tmpClus.addHit(h);
+                    }
+                }
+                muonmip.add(tmpClus); 
+            }    
+        }else{
+            List<Cluster> output = new Vector<Cluster>(); 
+            // Look for structure
+            List<Cluster> mipsOldInThisTree = findMipsOld(hitsInTree);
+            List<Cluster> mipsNewInThisTree = findMipsNew(hitsInTree);
+            List<Cluster> clumpsInThisTree = findClumps(hitsInTree);
+            List<Cluster> tightNeighbourClustersInThisTree = findNNClusters(hitsInTree);
+            // Treat the NN clusters as clumps for now
+            clumpsInThisTree.addAll(tightNeighbourClustersInThisTree);
+            // Record what we found
+            output.addAll(mipsOldInThisTree);
+            output.addAll(mipsNewInThisTree);
+            output.addAll(clumpsInThisTree);
+
+            if (debug) { System.out.println("size of output= " + output.size()); }
+            // Simple way by looking at direction of two clusters
+            List<Cluster> removed = new Vector<Cluster>();
+            if(output.size() > 1){
+                for(int i = 0 ; i < output.size() ; i ++){
+                    Cluster cluster1 = output.get(i);
+                    if(removed.contains(cluster1)) continue;
+                    Hep3Vector v1 = getClusterDirection(output.get(i));
+                    BasicCluster tmpClus = new BasicCluster();
+                    for(CalorimeterHit hit : cluster1.getCalorimeterHits()){
+                       tmpClus.addHit(hit);
+                    }
+
+                    for(int j = 1 ; j < output.size()-1 ; j ++){
+                        Cluster cluster2 = output.get(j);
+                        if(removed.contains(cluster2)) continue;
+                        Hep3Vector v2 = getClusterDirection(output.get(j));
+                        double cos = -1;
+                        if( v1 != null && v2 !=null){ cos =  VecOp.dot(v1, v2); }
+                        if(cos>0.8) {
+                            for(CalorimeterHit hit : cluster2.getCalorimeterHits()){
+                            tmpClus.addHit(hit);
+                            }
+                            removed.add(cluster2);
+                        } 
+                    }
+                    muonmip.add(tmpClus);
+                    removed.add(cluster1);
+                }    
+            }else { muonmip.addAll(output); }
+        }
+
+        if(debug){
+            int hitsize = 0 ;
+            for(Cluster clus : muonmip){
+                hitsize = hitsize + clus.getCalorimeterHits().size();
+                for(CalorimeterHit hit : clus.getCalorimeterHits()){
+                    examPosition(hit);
+                }
+            }
+        }
+
+        return muonmip;
+    }
+
+
+    public List<Cluster> createClusters(Set<CalorimeterHit> input_allhits)
+    {
+        List<Cluster> clusters = new Vector<Cluster>();
+        Set<CalorimeterHit> allhits = input_allhits;
+        for(CalorimeterHit hit : allhits){
+            BasicCluster tmpClus = new BasicCluster();
+            tmpClus.addHit(hit);
+            clusters.add(tmpClus);
+        }
+        return clusters;
+    }
+
+   protected void initTrackMatch() {
+        // Track-matching is complex. Use several matchers...
+        // Try matching with local helix extrap to MIP or generic cluster:
+        LocalHelixExtrapolationTrackMIPClusterMatcher mipMatch = new LocalHelixExtrapolationTrackMIPClusterMatcher(m_extrap);
+        LocalHelixExtrapolationTrackClusterMatcher genMatch = new LocalHelixExtrapolationTrackClusterMatcher(m_extrap);
+        DualActionTrackClusterMatcher localHelixMatchers = new DualActionTrackClusterMatcher(mipMatch, genMatch);
+        add(mipMatch);
+        add(genMatch);
+        m_trackClusterMatcher = localHelixMatchers;
+    }
+
+     protected List<Cluster> findMipsOld(HitMap unusedHits) {
+        Clusterer oldMipfinder = new TrackClusterDriver();
+        List<Cluster> mipClustersOld = oldMipfinder.createClusters(unusedHits);
+        //for (Cluster mip : mipClustersOld) {
+        //    for (CalorimeterHit hit : mip.getCalorimeterHits()) {
+        //        hitsInTree.remove(hit.getCellID());
+        //    }
+        //}
+        return mipClustersOld;
+    }
+                                                                                                                              
+    protected List<Cluster> findMipsNew(HitMap unusedHits) {
+        Clusterer newMipFinder = new org.lcsim.recon.cluster.mipfinder.NonProjectiveMipFinder(m_newMipFinderRadius);
+        List<Cluster> mipClustersNew = newMipFinder.createClusters(unusedHits);
+        //for (Cluster mip : mipClustersNew) {
+        //    for (CalorimeterHit hit : mip.getCalorimeterHits()) {
+        //        hitsInTree.remove(hit.getCellID());
+        //    }
+        //}
+        return mipClustersNew;
+    }
+                                                                                                                              
+    protected List<Cluster> findClumps(HitMap unusedHits) {
+        Clusterer clumpfinder = new ClumpFinder();
+        List<Cluster> clumpClusters = clumpfinder.createClusters(unusedHits);
+        for (Cluster clump : clumpClusters) {
+            if (clump.getCalorimeterHits().size() == 0) { throw new AssertionError("clump has no hits"); }
+            if (clump.getCalorimeterHits().contains(null)) { throw new AssertionError("null hit in clump"); }
+         //   for (CalorimeterHit hit : clump.getCalorimeterHits()) {
+         //       hitsInTree.remove(hit.getCellID());
+         //   }
+        }
+        return clumpClusters;
+    }
+
+     protected List<Cluster> findNNClusters(HitMap unusedHits) {
+        if (m_dU>0 && m_dV>0 && m_dL>0) {
+            int min_cells = m_minHitsToBeTreatedAsCluster/2;
+            Clusterer neighbourClusterer = new org.lcsim.recon.cluster.nn.NearestNeighborClusterer(m_dU, m_dV, m_dL, min_cells, 0.0);
+            List<Cluster> output = neighbourClusterer.createClusters(unusedHits);
+            for (Cluster clus : output) {
+                if (clus.getCalorimeterHits().size() == 0) { throw new AssertionError("clump has no hits"); }
+                if (clus.getCalorimeterHits().contains(null)) { throw new AssertionError("null hit in clump"); }
+                //for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+                //    hitsInTree.remove(hit.getCellID());
+                //}
+            }
+            return output;
+        } else {
+            // Do nothing -- return empty vector
+            return(new Vector<Cluster>());
+        }
+    }
+
+    protected void removePoorQualityMips(Collection<Cluster> mipList) {
+        SimpleMipQualityDecision check = new SimpleMipQualityDecision();
+        List<Cluster> mipsToRemove = new Vector<Cluster>();
+        for (Cluster mip : mipList) {
+            boolean passesCheck = check.valid(mip);
+            if (!passesCheck) {
+                mipsToRemove.add(mip);
+                if (debug) { System.out.println("Removed due to poor quality"); }
+            }
+        }
+        mipList.removeAll(mipsToRemove);
+    }
+
+                                                                                                                              
+    //examine neighbour hits
+    protected void examPosition(CalorimeterHit hit){
+	Set<Long> allHitIDs = new HashSet<Long>();
+	for (CalorimeterHit allhit : m_allhits)  allHitIDs.add(allhit.getCellID());
+
+	IDDecoder id = hit.getIDDecoder();
+	id.setID(hit.getCellID());
+	int count[] = new int[4];
+	for (int i =0 ; i < 4 ; i++){
+	    long[] neighbours = id.getNeighbourIDs(0,i+1,i+1);
+	    for(long neighbourID : neighbours)
+		if (allHitIDs.contains(neighbourID)) count[i]++;
+	}
+	if ( count[2] > 1 ) exam++ ; // if number of neighbour hits in 5x5 is larger than 1, count 1 
+	else exam = 0; // if number of neighbour is none, set to zero and start again from 0 
+
+	if(debug){
+	    Hep3Vector cellv = hitCellPosition(hit);
+            Hep3Vector v = new BasicHep3Vector(hit.getPosition());
+	    int p = (int) (v.magnitude());
+            int r = (int) Math.sqrt(v.x()*v.x() + v.y()*v.y());
+            int z = (int) v.z();
+            int cellp = (int) (cellv.magnitude());
+            int cellr = (int) Math.sqrt(cellv.x()*cellv.x() + cellv.y()*cellv.y());
+            int cellz = (int) cellv.z();
+	    System.out.println(" hit : p= " + p +" r=" + r + " , z " + z + " " + count[0] + " hits in 3x3, " + count[1] + " at " + id.getLayer() + " layer of " + hit.getSubdetector().getName());
+            System.out.println(" cell: p= " + cellp +" r=" + cellr + " , z " + cellz);
+	}
+        
+    }
+
+    protected  void removeUpperSensitiveLayer(Set<CalorimeterHit>  muhits){
+       Set<CalorimeterHit> upperhits = new HashSet<CalorimeterHit>();
+       for(CalorimeterHit hit : muhits){
+           String subdetName = hit.getSubdetector().getName();
+           Hep3Vector cellv = hitCellPosition(hit);
+           Hep3Vector v = new BasicHep3Vector(hit.getPosition());
+           int p = (int) (v.magnitude());
+           int r = (int) Math.sqrt(v.x()*v.x() + v.y()*v.y());
+           int z = (int) v.z();
+           int cellp = (int) (cellv.magnitude());
+           int cellr = (int) Math.sqrt(cellv.x()*cellv.x() + cellv.y()*cellv.y());
+           int cellz = (int) cellv.z();
+           boolean upper = true;
+           if(subdetName.contains("MuonEndcap")) upper = Math.abs(z) > Math.abs(cellz);  
+           if(subdetName.contains("MuonBarrel")) upper = r > cellr;
+           if(upper) upperhits.add(hit);
+       }
+       muhits.removeAll(upperhits);
+    }
+
+    //give the vector of hit
+    protected Hep3Vector hitCellPosition(CalorimeterHit hit){
+	IDDecoder id = hit.getIDDecoder();
+	id.setID(hit.getCellID());
+	double x = id.getPosition()[0];
+	double y = id.getPosition()[1];
+	double z = id.getPosition()[2];
+	Hep3Vector v = new BasicHep3Vector(x,y,z);
+	return v;
+    }
+
+    protected int getLayer(CalorimeterHit hit) {
+        org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+        id.setID(hit.getCellID());
+        int layer = id.getLayer();
+        return layer;
+    }
+
+    protected Hep3Vector getClusterDirection(Cluster clus){
+        Hep3Vector trackMuonUnit = new BasicHep3Vector();
+        if( clus.getCalorimeterHits().size() > 1){
+            CalorimeterHit muonhit1 = clus.getCalorimeterHits().get(clus.getCalorimeterHits().size()-1);
+            CalorimeterHit muonhit0 = clus.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);
+            }
+            trackMuonUnit = VecOp.unit(trackMuon);
+        } else { trackMuonUnit = null;}
+
+        return trackMuonUnit;
+    }
+}

lcsim/src/org/lcsim/recon/cluster/muonfinder
MuonFinder.java added at 1.1
diff -N MuonFinder.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MuonFinder.java	22 Oct 2008 23:24:34 -0000	1.1
@@ -0,0 +1,358 @@
+package org.lcsim.recon.cluster.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;
+
+/**
+ * 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.
+ *
+ * @author [log in to unmask]
+ * @version $Id: MuonFinder.java,v 1.1 2008/10/22 23:24:34 mcharles Exp $
+ */
+
+public class MuonFinder extends Driver{
+    protected String _tracklist;
+    protected String _outHitMap;
+    protected String _outMapTrackToClusters; 
+    protected String _inMu;
+    protected String _inCal;
+    protected double _bestmatch = 0.8;
+    protected boolean _debug = false;
+    protected boolean _useFirstTwoLayer = true;
+    EventHeader m_event;
+
+    protected HelixExtrapolator m_findCluster = new org.lcsim.recon.pfa.identifier.LocalHelixExtrapolator();
+    protected boolean useMuonHitsForDebug = false;
+
+    public MuonFinder(String tracklist, String inCal, String inMu, String outMap, String outMapTrackToClusters){
+        _inCal = inCal;
+        _inMu = inMu;
+        _outHitMap = outMap;
+        _tracklist = tracklist;
+	_outMapTrackToClusters = outMapTrackToClusters;
+        org.lcsim.contrib.uiowa.CalorimeterHitTimeCutDecision dec = new org.lcsim.contrib.uiowa.CalorimeterHitTimeCutDecision(100);
+        add(new ListFilterDriver(dec, "MuonBarrHits", "CorrMuonBarrHits", CalorimeterHit.class));
+        add(new ListFilterDriver(dec, "MuonEndcapHits", "CorrMuonEndcapHits", CalorimeterHit.class));
+        add(new SubsetFlagDriver("CorrMuonBarrHits"));
+        add(new SubsetFlagDriver("CorrMuonEndcapHits"));
+    }
+
+    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(useMuonHitsForDebug){
+            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, "CorrMuonBarrHits");
+            List<CalorimeterHit> allHitsMUEndcap = event.get(CalorimeterHit.class, "CorrMuonEndcapHits");
+            mudethits.addAll(allHitsMUBarrel);
+            mudethits.addAll(allHitsMUEndcap);
+        }else{
+            calhits.addAll(calHitMap.values());
+            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);
+
+        //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(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>>();
+        if(_debug){
+            System.out.println("Muon mip size= " + muonmips.size()+" (should be 1 for signle sample)");
+        }
+        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 = m_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){
+                                tpoint = tBarrel;
+                            }else{
+                                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.
+			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){
+                            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); 
+        event.put(_outHitMap, muonMap);
+        event.put(_outMapTrackToClusters, MuonMapTrackToCluster);
+        if(_debug){
+            List<Track> muonTracks = new Vector<Track>(MuonMapTrackToCluster.keySet());
+            List<Cluster> muonClusters = new Vector<Cluster>(MuonMapTrackToCluster.values());
+            int flag = 1<<LCIOConstants.CLBIT_HITS;
+            event.put("MuonClusters", muonClusters, Cluster.class, flag );
+            event.put("MuonList", muonTracks);
+        }
+    }
+
+    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();
+    }
+ 
+}

lcsim/src/org/lcsim/recon/cluster/muonfinder
MuonFinderWrapper.java added at 1.1
diff -N MuonFinderWrapper.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MuonFinderWrapper.java	22 Oct 2008 23:24:34 -0000	1.1
@@ -0,0 +1,109 @@
+package org.lcsim.recon.cluster.muonfinder;
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.event.util.*;
+import org.lcsim.util.*;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.HitInMCALDecision;
+import org.lcsim.util.hitmap.HitMap;
+
+/**
+ * Wrapper class for MuonFinder.
+ *
+ * @author [log in to unmask]
+ * @version $Id: MuonFinderWrapper.java,v 1.1 2008/10/22 23:24:34 mcharles Exp $
+ */
+
+public class MuonFinderWrapper extends Driver {
+
+    // Inputs
+    protected String m_trackList;
+    protected String m_inHitMap;
+    
+    // Outputs
+    protected String m_outMuonTrackClusterMap;
+    protected String m_outHitMap;
+    protected String m_outTrackList;
+
+    // Old/temporary
+    protected String _outHitMap;
+    protected String _outMapTrackToClusters;
+    protected String _inMu;
+    protected String _inCal;
+
+    // Option to skip
+    protected boolean m_debugSkip = false;
+    public void skip() { m_debugSkip = true; }
+
+    public MuonFinderWrapper(String inTrackList, String inHitMap, String outMuonTrackClusterMap,  String outHitMap, String outTrackList) 
+    {
+	// Inputs/outputs
+	m_trackList = inTrackList;
+	m_inHitMap = inHitMap;
+	m_outMuonTrackClusterMap = outMuonTrackClusterMap;
+	m_outHitMap = outHitMap;
+	m_outTrackList = outTrackList;
+	// Temporary stuff for old process routine
+	_inCal = new String("tmpCal");
+	_inMu = new String("tmpMu");
+	_outHitMap = new String("tmpMap");
+	_outMapTrackToClusters = new String("tmpMap2");
+    }
+
+    public void process(EventHeader event) {
+	// Read in
+	HitMap inHitMap = (HitMap)(event.get(m_inHitMap));
+	List<Track> inTrackList = event.get(Track.class, m_trackList);
+
+	// Massage inputs into the form wanted by old process routine
+	HitMap muonHits = new HitMap();
+	HitMap otherHits = new HitMap();
+	HitInMCALDecision dec = new HitInMCALDecision();
+	for (CalorimeterHit hit : inHitMap.values()) {
+	    long id = hit.getCellID();
+	    if (dec.valid(hit)) {
+		muonHits.put(id, hit);
+	    } else {
+		otherHits.put(id, hit);
+	    }
+	}
+	event.put(_inMu, muonHits);
+	event.put(_inCal, otherHits);
+
+	// Prepare outputs
+	List<Track> outTrackList = new Vector<Track>();
+	Map<Track,Cluster> outMuonTrackClusterMap = new HashMap<Track,Cluster>();
+	HitMap outHitMap = new HitMap(inHitMap);
+
+	if (m_debugSkip) {
+	    // Optionally, skip the muon finding
+	    outHitMap = inHitMap;
+	    outTrackList = inTrackList;
+	} else {
+	    // Call old process routine to do the work
+	    MuonFinder tmpWrappedMuonFinder = new MuonFinder(m_trackList, _inCal, _inMu, _outHitMap, _outMapTrackToClusters);
+	    tmpWrappedMuonFinder.process(event);
+	    outMuonTrackClusterMap = (Map<Track,Cluster>)(event.get(_outMapTrackToClusters));
+
+	    // Make HitMap of remaining hits:
+	    for (Cluster muonClus : outMuonTrackClusterMap.values()) {
+		for (CalorimeterHit hit : muonClus.getCalorimeterHits()) {
+		    long id = hit.getCellID();
+		    outHitMap.remove(id);
+		}
+	    }
+
+	    // Make List of remaining tracks:
+	    outTrackList.addAll(inTrackList);
+	    outTrackList.removeAll(outMuonTrackClusterMap.keySet());
+	}
+
+
+	// Write out
+	event.put(m_outMuonTrackClusterMap, outMuonTrackClusterMap);
+	event.put(m_outHitMap, outHitMap);
+	event.put(m_outTrackList, outTrackList);
+    }
+}
+

lcsim/src/org/lcsim/recon/cluster/muonfinder
SimpleMipQualityDecision.java added at 1.1
diff -N SimpleMipQualityDecision.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ SimpleMipQualityDecision.java	22 Oct 2008 23:24:34 -0000	1.1
@@ -0,0 +1,78 @@
+package org.lcsim.recon.cluster.muonfinder;
+
+import java.util.*; 
+
+import org.lcsim.util.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.*;
+import org.lcsim.event.*;
+
+/**
+ * A class to look at whether a calorimeter cluster is
+ * consistent with being a track/MIP. Criteria are based on
+ *
+ * modified from the code written by Mat.
+ *
+ * @author [log in to unmask]
+ * @version $Id: SimpleMipQualityDecision.java,v 1.1 2008/10/22 23:24:34 mcharles Exp $
+ */
+ 
+
+public class SimpleMipQualityDecision implements DecisionMakerSingle<Cluster>
+{
+    private int count;
+    boolean removeUpperHit = true;
+    public SimpleMipQualityDecision() {
+	// No user-settable options for now
+    }
+
+    public boolean valid(Cluster mip) {
+	// How many layers in this MIP?
+	// And how many instances of >1 hit in a layer?
+	Set<Integer> layersSeen = new HashSet<Integer>();
+        Set<Integer> isolayer = new HashSet<Integer>();
+        Map<Integer, Set<CalorimeterHit>> exam = new HashMap<Integer, Set<CalorimeterHit>>();
+	for (CalorimeterHit hit : mip.getCalorimeterHits()) {
+            String subdetName = hit.getSubdetector().getName();
+	    int layer = getLayer(hit);
+            //To avoid overlap with the other detector's layer
+            if(subdetName.contains("MuonEndcap")) layer = layer+100; 
+	    Integer layerInt = new Integer(layer);
+            if (!layersSeen.contains(layerInt)) {
+                layersSeen.add(layerInt);
+                isolayer.add(layerInt);
+                Set<CalorimeterHit> hits = exam.get(layerInt);
+                if(hits == null){
+                    hits = new HashSet<CalorimeterHit>();
+                    exam.put(layerInt,hits);
+                }
+                hits.add(hit);
+            } else {
+                Set<CalorimeterHit> hits = exam.get(layerInt);
+                if(hits == null){
+                    hits = new HashSet<CalorimeterHit>();
+                    exam.put(layerInt,hits);
+                }
+                hits.add(hit);
+                if(removeUpperHit){
+                    if(exam.get(layerInt).size() > 1){  isolayer.remove(layerInt); }
+                }else{
+                    if(exam.get(layerInt).size() > 2){  isolayer.remove(layerInt); }
+                }
+            }
+	}
+       
+        int count = 3;
+
+        boolean pass = (isolayer.size() >= count);
+        return pass;
+    }
+
+    protected int getLayer(CalorimeterHit hit) {
+        org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+        id.setID(hit.getCellID());
+        int layer = id.getLayer();
+        return layer;
+    }
+}
CVSspam 0.2.8