diff -N ShowerPointFinder.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ShowerPointFinder.java	13 Jul 2008 23:31:31 -0000	1.1
@@ -0,0 +1,268 @@

+package org.lcsim.recon.cluster.mipfinder;
+
+import java.util.*; 
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import hep.physics.vec.*;
+
+import java.io.IOException;
+import hep.physics.particle.properties.*;
+import org.lcsim.util.*;
+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 hep.physics.particle.Particle;
+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.util.swim.Line;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.*;
+import org.lcsim.util.decision.*;
+
+    public class ShowerPointFinder{
+    protected LocalHelixExtrapolator m_findCluster;
+    protected Map<Track, Set<Cluster>> m_newMapTrackToShowerComponents;
+    protected List<Cluster> m_mips;
+    protected Set<CalorimeterHit> m_allhits;
+    protected int exam = 0; //to test neighbour
+    protected boolean debug = false;
+
+    //sid01
+    protected double ecalRs=1270; //ECAL inner R
+    protected double ecalRe=1410; //ECAL outer R
+    protected double hcalRs=1410; //HCAL inner R
+    protected double hcalRe=2362; //HCAL outer R
+    protected double ecalZs=1680; //ECAL inner Z
+    protected double ecalZe=1820; //ECAL outer Z
+    protected double hcalZs=1820; //HCAL inner Z
+    protected double hcalZe=2772; //HCAL inner Z 
+
+    public ShowerPointFinder(LocalHelixExtrapolator findCluster, Map<Track, Set<Cluster>> newMapTrackToShowerComponents, List<Cluster> mips, Set<CalorimeterHit> allHits) {
+		m_findCluster = findCluster;
+		m_newMapTrackToShowerComponents = newMapTrackToShowerComponents;
+		m_mips = mips;
+		m_allhits = allHits; 
+    }
+
+	 
+    public Map<Track, BasicCluster> findMips(List<Track> tracksSortedByMomentum) {
+	Map<Track, BasicCluster> MapTrkToMIP = new HashMap<Track, BasicCluster>();
+	BasicCluster mipclus = new BasicCluster();
+	for(Track tr : tracksSortedByMomentum){ 
+	    Set<Cluster> clusoftrk = m_newMapTrackToShowerComponents.get(tr);
+	    Set<Cluster> mipsoftrk = new HashSet<Cluster>();
+	    for(Cluster cl : clusoftrk){
+		if(m_mips.contains(cl))	{
+		    mipsoftrk.add(cl);
+		}
+	    }
+	
+	    if(mipsoftrk.size() == 0) {//No any mip found
+		MapTrkToMIP.put(tr,mipclus);
+		continue; //go to next track
+	    }
+	    
+	    //new cluster for finding shower point
+	    Hep3Vector interceptPoint = m_findCluster.performExtrapolation(tr); //not need except debuging
+	    //Hep3Vector tangent =  m_findCluster.getTangent();
+	    //Hep3Vector tangentUnit = VecOp.unit(tangent);
+				
+	    if(debug){
+	    	Hep3Vector v = new BasicHep3Vector(tr.getMomentum());
+		System.out.println("Given track");
+		System.out.println("Debug: Size of Cluster of track= " + clusoftrk.size());
+		System.out.println("Debug: Size of Cluster of Mip of track= " + mipsoftrk.size());
+		System.out.println("Debug: trk momentum= " + v.magnitude());
+		if(interceptPoint != null) System.out.println("Debug: extra trk: pos= " + interceptPoint.magnitude() + " ("+ interceptPoint.x() + " " + interceptPoint.y() + " " + interceptPoint.z() + ")" );
+	    }
+
+	    //find closest mip to track
+	    SortedMap<Double, Cluster> mipsSortedByDistance = new TreeMap();
+	    for(Cluster cl : mipsoftrk){
+		CalorimeterHit firstHit =  cl.getCalorimeterHits().get(0);
+		Hep3Vector first = hitPosition(firstHit);
+		double distance = first.magnitude();
+		mipsSortedByDistance.put(distance, cl);
+	    }
+		
+	    //arrange by the length of the hit position
+	    SortedMap<Double, CalorimeterHit> sortedMIPbyPos = new TreeMap();
+	    Cluster cl = mipsSortedByDistance.get(mipsSortedByDistance.firstKey());
+    	    for( int i = 0 ; i < cl.getCalorimeterHits().size() ; i++){
+		CalorimeterHit hit=cl.getCalorimeterHits().get(i);
+		Hep3Vector v = hitPosition(hit);
+		double pmag = v.magnitude();
+            	sortedMIPbyPos.put(pmag, hit);
+	    }
+
+	    //test if the mip is enough closet to track otherwise take it as no mip.
+	    CalorimeterHit firstmip = sortedMIPbyPos.get(sortedMIPbyPos.firstKey());
+	    IDDecoder firstmipid = firstmip.getIDDecoder();
+	    firstmipid.setID(firstmip.getCellID());
+	    if(firstmipid.getLayer() > 3) { // if the mip is not in first 3 layers, put 0 in mipclus
+		MapTrkToMIP.put(tr,mipclus);
+		continue; //go to next track
+	    }
+	
+	    if(debug) System.out.println("Debug: Starting mip");
+	    Iterator itermip = sortedMIPbyPos.entrySet().iterator();
+	    for(int i=0 ; i < sortedMIPbyPos.size() ;i++){ 
+		Map.Entry entry = (Map.Entry) itermip.next();
+		Double key = (Double)entry.getKey();
+		CalorimeterHit hit = (CalorimeterHit)entry.getValue();
+		examPosition(hit);
+		mipclus.addHit(hit);
+	    }
+	    if(debug) System.out.println("Debug: The end of the mip" );
+	
+	    int lastIterationWithFoundHit = -1;
+	    int maxSkippedLayers = 3;
+	    for( int iIteration = 0 ; iIteration < 70 ; iIteration++){ // 70? until finding shower starting point
+		//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();
+
+		// 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, 2, 2);
+		for (long lastnearID : nearbyHitArray ) { nearIDs.add(lastnearID); }
+
+		Hep3Vector last0 = hitPosition(hit);
+		CalorimeterHit hit1 = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-2);
+		Hep3Vector last1 = hitPosition(hit1);
+
+		//Searh next hit in next layer and arrange in angular order
+		SortedMap<Double, CalorimeterHit> sortedHitbyPos= new TreeMap();
+		for (CalorimeterHit hitall : m_allhits){
+		    IDDecoder idall = hitall.getIDDecoder();
+		    idall.setID(hitall.getCellID());
+		    double xall = idall.getPosition()[0];
+		    double yall = idall.getPosition()[1];
+		    double zall = idall.getPosition()[2];
+		    double hitmag = Math.sqrt(xall*xall+yall*yall+zall*zall);
+		    Hep3Vector hitpos = new BasicHep3Vector(xall,yall,zall);
+		    int hitLayer = idall.getLayer();
+		    //use geometric information and layer number to find next hit
+		    boolean EcalToHcal= lastlayer == 30 && hitmag > p && idall.getLayer() == 0;
+		    boolean InSide = (hitmag > p || hitLayer > lastlayer) && nearIDs.contains(hitall.getCellID());
+		    boolean nearInnerEdgeBarrelEcal = (r < ecalRs && r > ecalRs-4);
+		    boolean nearOuterEdgeEndcapEcal = (Math.abs(z) < ecalZe && Math.abs(z) > ecalZs);
+		    boolean nearInnerEdgeBarrelHcal = (r < hcalRs && r > hcalRs-4);
+		    boolean nearOuterEdgeEndcapHcal = (Math.abs(z) < hcalZe && Math.abs(z) > hcalZs);		    
+		    boolean EndToBarrelEcal = nearInnerEdgeBarrelEcal && nearOuterEdgeEndcapEcal && (hitmag > p) && (idall.getLayer()==0);
+		    boolean EndToBarrelHcal = nearInnerEdgeBarrelHcal && nearOuterEdgeEndcapHcal && (hitmag > p) && (idall.getLayer()==0);
+
+		    if(debug){
+			if(EndToBarrelEcal || EndToBarrelHcal ) System.out.println("This is a boundary condition");
+		    }				
+
+		    if(EcalToHcal || InSide || EndToBarrelEcal || EndToBarrelHcal){
+			Hep3Vector test = VecOp.sub(hitpos, last0);
+			Hep3Vector test_unit = VecOp.unit(test);
+			Hep3Vector lastmip = VecOp.sub(last0, last1);
+			Hep3Vector lastmip_unit = VecOp.unit(lastmip);
+			double a = Math.acos(VecOp.dot(lastmip_unit, test_unit));
+			sortedHitbyPos.put(a, hitall);
+		    }		
+
+		}		
+		
+		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 > 2) {
+		    if(debug) System.out.println("Debug: showering is already started " + exam + " layers ago" );
+		    if(mipclus.getCalorimeterHits().size() < 5) break; // if hit size is fewer than 5, don't remove 
+		    else for(int k=0; k < 3 ; 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() < 4) k = mipclus.getCalorimeterHits().size();
+		else k = 4;
+		for(int i =0; i < k ;i++){
+		    CalorimeterHit hit = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-1-i);
+		    Hep3Vector v = hitPosition(hit);
+		    System.out.println("Debug: last "+i+" pos= "+v.magnitude()+" (" +v.x() + " " + v.y() + " " + v.z()+ ")" );
+		}
+		System.out.println("MIP hit size= " + mipclus.getCalorimeterHits().size());
+	    }		
+
+	    MapTrkToMIP.put(tr,mipclus); 
+	} //track loop
+
+	return MapTrkToMIP;
+    }
+
+    //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[3] >=2 ) exam++ ; // if number of neighbour hits in 9x9 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 v = hitPosition(hit);
+	    int p = (int) (v.magnitude());
+	    System.out.println(" hit: pos=" + p + " , " + count[0] + " hits in 3x3, " + count[1] + " hits in 5x5 " + count[2] + " hits in 7x7 " + count[3] + " hits in 9x9 " + id.getLayer() + " layer of " + hit.getSubdetector().getName());
+	}
+    }
+
+    //give the vector of hit
+    protected Hep3Vector hitPosition(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;
+    }
+}