diff -N MergeClustersCrossingSubDetectorBoundaries.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MergeClustersCrossingSubDetectorBoundaries.java	28 Sep 2010 10:16:35 -0000	1.1
@@ -0,0 +1,440 @@

+package org.lcsim.recon.pfa.structural;
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.util.*;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.geometry.*;
+import hep.physics.vec.*;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.util.CalorimeterInformation;
+import org.lcsim.geometry.Calorimeter.CalorimeterType;
+
+/**
+ * An algorithm to merge clusters that get split because crossing 
+ * subdetector boundaries.
+ *
+ * @version $Id: MergeClustersCrossingSubDetectorBoundaries.java,v 1.1 2010/09/28 10:16:35 zaidan Exp $
+ * @author Remi Zaidan <[log in to unmask]>
+ */
+
+public class MergeClustersCrossingSubDetectorBoundaries extends Driver {
+
+    protected String m_inputName;
+    protected String m_outputName;
+
+    protected List<Cluster> m_inputLargeClusterList;
+    protected List<Cluster> m_outputLargeClusterList;
+
+    protected double m_minOverlapAllowed = 0;
+    protected int m_minHitsPerCluster = 2;
+
+    protected EventHeader m_event;
+    protected boolean m_debug;
+    protected boolean m_init;
+
+    protected int m_EcalBarrelNLayers;
+    protected int m_EcalEndcapNLayers;
+
+    protected String m_EcalBarrelName;
+    protected String m_EcalEndcapName;
+    protected String m_HcalBarrelName;
+    protected String m_HcalEndcapName;
+    
+    protected double m_maxEcalEndcapR;
+    protected double m_maxHcalEndcapR;
+    protected double m_maxEcalBarrelZ;
+
+    protected double m_EcalEndcapROffset = 0;
+    protected double m_HcalEndcapROffset = 0;
+    protected double m_EcalBarrelZOffset = 0;
+    
+    public enum ConnectionType { 
+	ECAL_BARREL_TO_HCAL_BARREL,
+        ECAL_ENDCAP_TO_HCAL_ENDCAP,
+	ECAL_ENDCAP_TO_HCAL_BARREL,
+	ECAL_BARREL_TO_HCAL_ENDCAP,
+	ECAL_BARREL_TO_ECAL_ENDCAP,
+        HCAL_BARREL_TO_HCAL_ENDCAP
+    }
+
+    public MergeClustersCrossingSubDetectorBoundaries(String input, String output){
+	m_inputName = input;
+	m_outputName = output;
+
+	m_debug = false;
+
+	m_init = false;
+    }
+
+    /**
+     * Loops on clusters in input list and match them according to their
+     * overlap at the Hcal/Ecal boundary.
+     * Overlapping clusters are merged.
+     * Output consists of the merged clusters and the remaining non overlapping
+     * clusters from input list.
+     */
+    public void process(EventHeader event) {
+	
+	m_event = event;
+	
+	if(!m_init){
+            m_init = true;
+            CalorimeterInformation ci = CalorimeterInformation.instance();
+	    m_EcalBarrelNLayers = ci.getNLayers(CalorimeterType.EM_BARREL);
+	    m_EcalEndcapNLayers = ci.getNLayers(CalorimeterType.EM_ENDCAP);
+	    m_EcalBarrelName = ci.getName(CalorimeterType.EM_BARREL);
+	    m_EcalEndcapName = ci.getName(CalorimeterType.EM_ENDCAP);
+	    m_HcalBarrelName = ci.getName(CalorimeterType.HAD_BARREL);
+	    m_HcalEndcapName = ci.getName(CalorimeterType.HAD_ENDCAP);
+	    m_maxEcalEndcapR = ci.getRMax(CalorimeterType.EM_ENDCAP);
+	    m_maxHcalEndcapR = ci.getRMax(CalorimeterType.HAD_ENDCAP);
+	    m_maxEcalBarrelZ = ci.getRMax(CalorimeterType.EM_BARREL);
+	    m_maxEcalEndcapR -= m_EcalEndcapROffset;
+	    m_maxHcalEndcapR -= m_HcalEndcapROffset;
+	    m_maxEcalBarrelZ -= m_EcalBarrelZOffset;
+	    if(m_debug){
+		System.out.println("Calorimeter info: ");
+		System.out.println("  m_EcalBarrelNLayers = "+m_EcalBarrelNLayers);
+		System.out.println("  m_EcalEndcapNLayers = "+m_EcalEndcapNLayers);
+		System.out.println("  m_EcalBarrelName = "+m_EcalBarrelName);
+		System.out.println("  m_EcalEndcapName = "+m_EcalEndcapName);
+		System.out.println("  m_HcalBarrelName = "+m_HcalBarrelName);
+		System.out.println("  m_HcalEndcapName = "+m_HcalEndcapName);
+		System.out.println("  m_maxEcalEndcapR = "+m_maxEcalEndcapR);
+		System.out.println("  m_maxHcalEndcapR = "+m_maxHcalEndcapR);
+		System.out.println("  m_maxEcalBarrelZ = "+m_maxEcalBarrelZ);
+	    }
+	}
+	
+	// get input list and copy it's contents to local list
+	List<Cluster> inputLargeClusterOriginalList = m_event.get(Cluster.class, m_inputName);
+	m_inputLargeClusterList = new Vector<Cluster>();
+	m_inputLargeClusterList.addAll(inputLargeClusterOriginalList);
+	
+	// sort clusters by size: Larger clusters sould come first
+	Collections.sort(m_inputLargeClusterList, new ClusterSizeSort());
+	
+	// find connections
+	m_outputLargeClusterList = new Vector<Cluster>();
+	Map<Cluster, Map<ConnectionType, Cluster> > connectionMap = new HashMap<Cluster, Map<ConnectionType, Cluster> >();
+	for(int i=0; i<m_inputLargeClusterList.size(); i++){
+	    Cluster clus1 = m_inputLargeClusterList.get(i);
+	    if(m_debug){
+		System.out.println("searching for clusters to merge with cluster ["+i+"] with "+clus1.getCalorimeterHits().size()+" hits");
+	    }
+	    if(clus1.getCalorimeterHits().size() < m_minHitsPerCluster) break;
+	    for(ConnectionType c : ConnectionType.values()){
+		if(m_debug){
+		    System.out.println("looking for connections across "+c+" border");
+		}
+		Map<ConnectionType, Cluster> clus1Connections = connectionMap.get(clus1);
+		if(clus1Connections != null && clus1Connections.get(c) != null){
+		    if(m_debug){
+			System.out.println("Cluster was already connected across "+c+" border");
+		    }
+		    continue;
+		}
+		double maxOverlap = m_minOverlapAllowed;
+		Cluster clusmatch = null;
+		for(int j=i+1; j<m_inputLargeClusterList.size(); j++){
+		    Cluster clus2 = m_inputLargeClusterList.get(j);
+		    if(clus2.getCalorimeterHits().size() < m_minHitsPerCluster) break;
+		    Map<ConnectionType, Cluster> clus2Connections = connectionMap.get(clus2);
+		    if(clus2Connections != null && clus2Connections.get(c) != null) continue;
+		    double overlap = getOverlapFactor(clus1, clus2, c);
+		    if(overlap > maxOverlap){
+			maxOverlap = overlap;
+			clusmatch = clus2;
+		    }
+		}
+		if(clusmatch == null){
+		    if(m_debug){
+			System.out.println("no cluster match found across "+c+" border");
+		    }
+		    continue;
+		}
+
+		if(m_debug){
+		    System.out.println("Found overlap on "+c+" border between cluster with "+clus1.getCalorimeterHits().size()+" hits and cluster with "+clusmatch.getCalorimeterHits().size()+" hits: overlap factor is "+maxOverlap);
+		}
+
+		if(clus1Connections == null){
+		    clus1Connections = new HashMap<ConnectionType, Cluster>();
+		}
+		clus1Connections.put(c, clusmatch);
+		connectionMap.put(clus1, clus1Connections);
+		
+		Map<ConnectionType, Cluster> clusmatchConnections = connectionMap.get(clusmatch);
+		if(clusmatchConnections == null){
+		    clusmatchConnections = new HashMap<ConnectionType, Cluster>();
+		}
+		clusmatchConnections.put(c, clus1);
+		connectionMap.put(clusmatch, clusmatchConnections);
+	    }
+	}
+	
+	// merge clusters
+	Set<Cluster> mergedClusters = new HashSet<Cluster>();
+	for(Cluster clus : m_inputLargeClusterList){
+	    if(mergedClusters.contains(clus)) continue;
+	    Set<Cluster> clusmatches = new HashSet<Cluster>();
+	    getClustersToMergeWith(clus, connectionMap, clusmatches);
+	    BasicCluster mergedClus = new BasicCluster();
+	    for(Cluster clusmatch : clusmatches){
+		mergedClus.addCluster(clusmatch);
+		mergedClusters.add(clusmatch);
+	    }
+	    m_outputLargeClusterList.add(mergedClus);
+	}
+	
+	// Store new output
+	int flag = 1<<LCIOConstants.CLBIT_HITS;
+	m_event.put(m_outputName, m_outputLargeClusterList, Cluster.class, flag);
+	m_event.put(m_inputName+"Original", m_inputLargeClusterList, Cluster.class, flag);
+    }
+    
+    /** 
+     * Gets overlap factor defined as the size of the overlap region
+     * at a given border in (phi x theta) devided by the size 
+     * of the (phi x theta) region on the inner side.
+     * This supposes that each of the input clusters has hits on only one side
+     * of the border while the other one has only hits hits on the other side.
+     * To define sensible overlap regions, clusters must be of a certain
+     * size: non zero span in both phi and theta on both sides of Ecal/Hcal
+     * border are required.
+     */
+    protected double getOverlapFactor(Cluster clus1, Cluster clus2, ConnectionType connectionType){
+	
+	List<CalorimeterHit> hits_11 = null;
+	List<CalorimeterHit> hits_12 = null;
+	List<CalorimeterHit> hits_21 = null;
+	List<CalorimeterHit> hits_22 = null;
+	
+	switch(connectionType){
+	    
+	case ECAL_BARREL_TO_HCAL_BARREL:
+	    // get hits on Ecal/Hcal Barrel border: last Ecal layer and first Hcal layer
+	    hits_11 = getHitsOnLayer(clus1, m_EcalBarrelName, m_EcalBarrelNLayers-1);
+	    hits_12 = getHitsOnLayer(clus2, m_EcalBarrelName, m_EcalBarrelNLayers-1);
+	    hits_21 = getHitsOnLayer(clus1, m_HcalBarrelName, 0);
+	    hits_22 = getHitsOnLayer(clus2, m_HcalBarrelName, 0);
+	    break;
+	    
+	case ECAL_ENDCAP_TO_HCAL_ENDCAP:
+	    // get hits on Ecal/Hcal Endcap border: last Ecal layer and first Hcal layer
+	    hits_11 = getHitsOnLayer(clus1, m_EcalEndcapName, m_EcalEndcapNLayers-1);
+	    hits_12 = getHitsOnLayer(clus2, m_EcalEndcapName, m_EcalEndcapNLayers-1);
+	    hits_21 = getHitsOnLayer(clus1, m_HcalEndcapName, 0);
+	    hits_22 = getHitsOnLayer(clus2, m_HcalEndcapName, 0);
+	    break;
+	    
+	case ECAL_BARREL_TO_ECAL_ENDCAP:
+	    // get hits on Ecal Barrel/Endcap border: first Ecal barrel layer and max Ecal Endcap radius
+	    hits_11 = getHitsOnEdge(clus1, m_EcalEndcapName, m_maxEcalEndcapR);
+	    hits_12 = getHitsOnEdge(clus2, m_EcalEndcapName, m_maxEcalEndcapR);
+	    hits_21 = getHitsOnLayer(clus1, m_EcalBarrelName, 0);
+	    hits_22 = getHitsOnLayer(clus2, m_EcalBarrelName, 0);
+	    break;
+	    
+	case HCAL_BARREL_TO_HCAL_ENDCAP:
+	    // get hits on Hcal Barrel/Endcap border: first Hcal barrel layer and max Hcal Endcap radius
+	    hits_11 = getHitsOnEdge(clus1, m_HcalEndcapName, m_maxHcalEndcapR);
+	    hits_12 = getHitsOnEdge(clus2, m_HcalEndcapName, m_maxHcalEndcapR);
+	    hits_21 = getHitsOnLayer(clus1, m_HcalBarrelName, 0);
+	    hits_22 = getHitsOnLayer(clus2, m_HcalBarrelName, 0);
+	    break;
+	    
+	case ECAL_ENDCAP_TO_HCAL_BARREL:
+	    // not possible
+	    return 0;
+	    
+	case ECAL_BARREL_TO_HCAL_ENDCAP:
+	    // get hits on Ecal Barrel/Hcal Endcap border: max Ecal Barrel z and first Hcal Endcap layer
+	    hits_11 = getHitsOnEdge(clus1, m_EcalBarrelName, m_maxEcalBarrelZ);
+	    hits_12 = getHitsOnEdge(clus2, m_EcalBarrelName, m_maxEcalBarrelZ);
+	    hits_21 = getHitsOnLayer(clus1, m_HcalEndcapName, 0);
+	    hits_22 = getHitsOnLayer(clus2, m_HcalEndcapName, 0);
+	    break;
+	}    
+	return getOverlapFactor(hits_11, hits_12, hits_21, hits_22);   
+    }
+	
+    protected double getOverlapFactor(List<CalorimeterHit> hits_in1, List<CalorimeterHit> hits_in2, List<CalorimeterHit> hits_out1, List<CalorimeterHit> hits_out2){
+
+	if(hits_in1 == null) return 0;
+	if(hits_in2 == null) return 0;
+	if(hits_out1 == null) return 0;
+	if(hits_out2 == null) return 0;
+
+	int in1 = hits_in1.size();
+	int in2 = hits_in2.size();
+	int out1 = hits_out1.size();
+	int out2 = hits_out2.size();
+
+	// check if clusters have hits on opposite sides of the border
+	if(! ( (in1==0 && in2!=0 && out1!=0 && out2==0) ||
+	       (in1!=0 && in2==0 && out1==0 && out2!=0) ) ) return 0;
+
+	List<CalorimeterHit> hits_in;
+	List<CalorimeterHit> hits_out;
+	if(in1==0 && in2!=0 && out1!=0 && out2==0){
+	    hits_in = hits_in2;
+	    hits_out = hits_out1;
+	}else{
+	    hits_in = hits_in1;
+	    hits_out = hits_out2;
+	}
+
+	// compute clusters span in phi and eta
+	double phimin_in = 999;
+	double phimax_in = -999;
+	double thetamin_in = 999;
+	double thetamax_in = -999;
+	for(CalorimeterHit hit : hits_in){
+	    SpacePoint pos = new SpacePoint(new BasicHep3Vector(hit.getPosition()));
+	    double phi = pos.phi();
+	    double theta = pos.theta();
+	    if(phi < phimin_in) phimin_in = phi;
+	    if(phi > phimax_in) phimax_in = phi;
+	    if(theta < thetamin_in) thetamin_in = theta;
+	    if(theta > thetamax_in) thetamax_in = theta;
+	}
+	if(phimin_in == -999) return 0;
+	if(phimax_in == 999) return 0;
+	if(thetamin_in == -999) return 0;
+	if(thetamax_in == 999) return 0;
+
+	// check for too small region
+	if(!(phimax_in - phimin_in > 0)) return 0;
+	if(!(thetamax_in - thetamin_in > 0)) return 0;
+	
+	double phimin_out = 999;
+	double phimax_out = -999;
+	double thetamin_out = 999;
+	double thetamax_out = -999;
+	for(CalorimeterHit hit : hits_out){
+	    SpacePoint pos = new SpacePoint(new BasicHep3Vector(hit.getPosition()));
+	    double phi = pos.phi();
+	    double theta = pos.theta();
+	    if(phi < phimin_out) phimin_out = phi;
+	    if(phi > phimax_out) phimax_out = phi;
+	    if(theta < thetamin_out) thetamin_out = theta;
+	    if(theta > thetamax_out) thetamax_out = theta;
+	}
+	if(phimin_out == -999) return 0;
+	if(phimax_out == 999) return 0;
+	if(thetamin_out == -999) return 0;
+	if(thetamax_out == 999) return 0;
+	
+	// check for too small region
+	if(!(phimax_out - phimin_out > 0)) return 0;
+	if(!(thetamax_out - thetamin_out > 0)) return 0;
+	
+	// define overlap region
+	double phimin = phimin_in > phimin_out ? phimin_in : phimin_out;
+	double phimax = phimax_in < phimax_out ? phimax_in : phimax_out;
+	double thetamin = thetamin_in > thetamin_out ? thetamin_in : thetamin_out;
+	double thetamax = thetamax_in < thetamax_out ? thetamax_in : thetamax_out;
+	
+	// check for overlap
+	if(!(phimax - phimin > 0)) return 0;
+	if(!(thetamax - thetamin > 0)) return 0;
+	
+	// overlap region normalized to inner span
+	double overlap_phi = (phimax - phimin) / (phimax_in - phimin_in);
+	double overlap_theta = (thetamax - thetamin) / (thetamax_in - thetamin_in);
+	
+	return overlap_phi * overlap_theta;
+    }
+    
+    /**
+     * Get list of hits on a given layer and detector
+     */
+    protected List<CalorimeterHit> getHitsOnLayer(Cluster clus, String det, int layer){
+	
+	List<CalorimeterHit> hits = clus.getCalorimeterHits();
+	List<CalorimeterHit> hitsOnLayer = new Vector<CalorimeterHit>();
+	
+	int i=0;
+	for(CalorimeterHit hit : hits){
+            IDDecoder id = hit.getIDDecoder();
+            id.setID(hit.getCellID());
+	    if(!id.getSubdetector().getName().equals(det)){
+		return hitsOnLayer;
+	    }
+	    if(id.getLayer() == layer){
+		hitsOnLayer.add(hit);
+	    }
+	    i++;
+	}
+	return hitsOnLayer;
+    }
+    
+    /**
+     * Get list of hits on a given edge of a detector: based on maximum z/R for Barrel/Endcap
+     */
+    protected List<CalorimeterHit> getHitsOnEdge(Cluster clus, String det, double xmax){
+	
+	List<CalorimeterHit> hits = clus.getCalorimeterHits();
+	List<CalorimeterHit> hitsOnEdge = new Vector<CalorimeterHit>();
+	
+	for(CalorimeterHit hit : hits){
+            IDDecoder id = hit.getIDDecoder();
+            id.setID(hit.getCellID());
+	    if(!id.getSubdetector().getName().equals(det)){
+		return hitsOnEdge;
+	    }
+	    double x = 0;
+	    BasicHep3Vector pos = new BasicHep3Vector(hit.getPosition());
+	    if(det.equals(m_EcalBarrelName) || det.equals(m_HcalBarrelName)){
+		x = Math.abs(pos.z());
+	    }
+	    if(det.equals(m_EcalEndcapName) || det.equals(m_HcalEndcapName)){
+		x = Math.sqrt(pos.x()*pos.x() +pos.y()*pos.y());
+	    }
+	    if(x >= xmax){
+		hitsOnEdge.add(hit);
+	    }
+	}
+	return hitsOnEdge;
+    }
+    
+    /**
+     * recursively find all clusters to be merged with a given cluster
+     */
+    protected void getClustersToMergeWith(Cluster clus, Map<Cluster, Map<ConnectionType, Cluster> > connectionMap, Set<Cluster> clusSet){
+	
+	if(clusSet.contains(clus)) return;
+	clusSet.add(clus);
+	
+	Map<ConnectionType, Cluster> connections = connectionMap.get(clus);
+
+	if(connections == null) return;
+	
+	for(ConnectionType c : connections.keySet()){
+	    Cluster clusmatch = connections.get(c);
+	    if(!clusSet.contains(clusmatch)){
+		getClustersToMergeWith(clusmatch, connectionMap, clusSet);
+	    }
+	}
+    }
+    
+    /**
+     * Larger clusters sould come first
+     */
+    protected class ClusterSizeSort implements Comparator<Cluster> {
+        public ClusterSizeSort() {}
+        public int compare(Cluster c1, Cluster c2) {
+	    int s1 = c1.getCalorimeterHits().size();
+	    int s2 = c2.getCalorimeterHits().size();
+            if (s1 < s2) {
+                return 1;
+            } else if (s1 > s2) {
+                return -1;
+            } else {
+                return 0;
+            }
+        }
+    }
+}