lcsim/src/org/lcsim/contrib/uiowa
diff -N MIPReassignmentAlgorithm.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MIPReassignmentAlgorithm.java 2 Jul 2008 00:07:42 -0000 1.1
@@ -0,0 +1,311 @@
+package org.lcsim.contrib.uiowa;
+
+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 MIPReassignmentAlgorithm implements ReassignClustersAlgorithm {
+ protected double m_limit;
+ protected LocalHelixExtrapolator m_findCluster;
+ protected Map<Track, Set<Cluster>> m_newMapTrackToShowerComponents;
+ protected List<Cluster> m_mips;
+ protected Set<CalorimeterHit> m_allhits;
+ static int exam = 0;
+ private boolean debug = false;
+
+ public MIPReassignmentAlgorithm(double limit, LocalHelixExtrapolator findCluster, Map<Track, Set<Cluster>> newMapTrackToShowerComponents, List<Cluster> mips, Set<CalorimeterHit> allHits) {
+ m_limit = limit;
+ m_findCluster = findCluster;
+ m_newMapTrackToShowerComponents = newMapTrackToShowerComponents;
+ m_mips = mips;
+ m_allhits = allHits;
+ }
+
+
+ public Double computeFigureOfMerit(Track tr, Cluster clus) {
+
+ List<Hep3Vector> lastTwoPositions = computePointsOfShower(tr);
+
+ Hep3Vector last0pos = lastTwoPositions.get(0);
+ Hep3Vector last1pos = lastTwoPositions.get(1);
+ Hep3Vector tangent = VecOp.sub(last0pos, last1pos);
+ Hep3Vector tangentUnit = VecOp.unit(tangent);
+ Hep3Vector clusterPosition = new BasicHep3Vector(clus.getPosition());
+ Hep3Vector displacement = VecOp.sub(clusterPosition, last1pos);
+ Hep3Vector displacementUnit = VecOp.unit(displacement);
+
+ double angle=Math.acos(VecOp.dot(tangentUnit, displacementUnit));
+ if (angle < m_limit) return new Double(angle);
+
+ return null;
+ }
+
+
+ protected List<Hep3Vector> computePointsOfShower(Track tr) {
+ 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(debug) System.out.println("Given track");
+
+ Hep3Vector interceptPoint = m_findCluster.performExtrapolation(tr);
+ Hep3Vector first = new BasicHep3Vector();
+ Hep3Vector trackThreeMomentum = new BasicHep3Vector(tr.getMomentum());
+ double trackMomentum = trackThreeMomentum.magnitude();
+
+ if(debug){
+ 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= " + trackMomentum);
+ if(interceptPoint != null) System.out.println("Debug: extra trk: pos= " + interceptPoint.magnitude() + " ("+ interceptPoint.x() + " " + interceptPoint.y() + " " + interceptPoint.z() + ")" );
+ }
+
+ SortedMap<Double, Cluster> mipsSortedByDistance = new TreeMap();
+ for(Cluster cl : mipsoftrk){
+ CalorimeterHit firstHit = cl.getCalorimeterHits().get(0);
+ IDDecoder id = firstHit.getIDDecoder();
+ id.setID(firstHit.getCellID());
+ double x = id.getPosition()[0];
+ double y = id.getPosition()[1];
+ double z = id.getPosition()[2];
+ first = new BasicHep3Vector(x, y, z);
+ Hep3Vector displacement = new BasicHep3Vector();
+ if (interceptPoint != null) {
+ displacement = VecOp.sub(first, interceptPoint);
+ double distance = displacement.magnitude();
+ mipsSortedByDistance.put(distance, cl);
+ }
+ }
+
+ SortedMap<Integer, CalorimeterHit> sortedMipBylayer = new TreeMap();
+ int t=0;
+ Iterator iter = mipsSortedByDistance.entrySet().iterator();
+ while(iter.hasNext() && t == 0){
+ t++;
+ Map.Entry entry = (Map.Entry) iter.next();
+ Double key = (Double)entry.getKey();
+ Cluster cl = (Cluster)entry.getValue();
+ for( int i = 0 ; i < cl.getCalorimeterHits().size() ; i++){
+ CalorimeterHit hit=cl.getCalorimeterHits().get(i);
+ IDDecoder idMIP = hit.getIDDecoder();
+ idMIP.setID(hit.getCellID());
+ int layerMIP = idMIP.getLayer();
+ sortedMipBylayer.put(layerMIP, hit);
+ }
+ }
+
+
+ Hep3Vector last0pos = new BasicHep3Vector();
+ Hep3Vector last1pos = new BasicHep3Vector();
+ Hep3Vector last2pos = new BasicHep3Vector();
+ Hep3Vector last3pos = new BasicHep3Vector();
+ Hep3Vector last4pos = new BasicHep3Vector();
+
+
+ Set<Long> allHitIDs = new HashSet<Long>();
+ for (CalorimeterHit hit : m_allhits) allHitIDs.add(hit.getCellID());
+
+ if(debug) System.out.println("Debug: Starting mip");
+
+ int lastlayer = 0;
+ long[] lastlayer5x5x2 = new long[26];
+ String lastDet = new String();
+ int sizeofmip = sortedMipBylayer.size();
+ Iterator itermip = sortedMipBylayer.entrySet().iterator();
+ for(int i=0 ; i < sizeofmip ;i++){
+ Map.Entry entry = (Map.Entry) itermip.next();
+ Integer key = (Integer)entry.getKey();
+ CalorimeterHit hit = (CalorimeterHit)entry.getValue();
+
+ Hep3Vector last= lastPosition(hit,i,allHitIDs);
+
+ if(i == (sizeofmip-1)) {
+ IDDecoder lastid = hit.getIDDecoder();
+ lastid.setID(hit.getCellID());
+ lastlayer = lastid.getLayer();
+ last0pos = last;
+ lastDet = hit.getSubdetector().getName();
+ lastlayer5x5x2 = lastid.getNeighbourIDs(1, 2, 2);
+ }
+ if(i == (sizeofmip-2)) last1pos = last;
+ if(i == (sizeofmip-3)) last2pos = last;
+ if(i == (sizeofmip-4)) last3pos = last;
+ }
+ int s=sizeofmip;
+ if(debug) System.out.println("Debug: The end of the mip" );
+
+
+ Hep3Vector lastv_mip = VecOp.sub(last0pos, last1pos);
+ Hep3Vector lastv_mip_unit = VecOp.unit(lastv_mip);
+
+ int cl = 0 ;
+ for( int i = lastlayer+1 ; i < lastlayer+100 ; i++){
+ SortedMap<Double, CalorimeterHit> sortedHitbyAngle= new TreeMap();
+ for (CalorimeterHit hit : m_allhits) {
+
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ double x = id.getPosition()[0];
+ double y = id.getPosition()[1];
+ double z = id.getPosition()[2];
+
+ Set<Long> nearIDs = new HashSet<Long>();
+ for (long lastnearID : lastlayer5x5x2 ) nearIDs.add(lastnearID);
+ boolean cond = ((id.getLayer() != 0) && (lastDet == hit.getSubdetector().getName()) && nearIDs.contains(hit.getCellID())) || ((id.getLayer() == 0) && (lastDet != hit.getSubdetector().getName()));
+
+ if(i <= 30) cl = i;
+ if(i > 30) cl = i-31;
+
+ if(cl==id.getLayer() && cond){
+ Hep3Vector hitpos = new BasicHep3Vector();
+ hitpos = new BasicHep3Vector(x,y,z);
+ Hep3Vector test = VecOp.sub(hitpos, last0pos);
+ double d = test.magnitude();
+ Hep3Vector test_unit = VecOp.unit(test);
+ double a = Math.acos(VecOp.dot(lastv_mip_unit, test_unit));
+ if ( a < 0.79) { // 45 degree = 0.79 radian
+ sortedHitbyAngle.put(a, hit);
+ IDDecoder idt = hit.getIDDecoder();
+ idt.setID(hit.getCellID());
+ }
+ }
+
+ }
+
+ Hep3Vector last = new BasicHep3Vector();
+ Iterator it = sortedHitbyAngle.entrySet().iterator();
+ int te = 0;
+ while(it.hasNext() && te == 0){
+ te++;
+ Map.Entry e = (Map.Entry) it.next();
+ Double k = (Double)e.getKey();
+ CalorimeterHit hit = (CalorimeterHit)e.getValue();
+ last=lastPosition(hit,s,allHitIDs);
+ s++;
+
+
+
+ if(!(last3pos == null)) last4pos = last3pos;
+ if(!(last2pos == null)) last3pos = last2pos;
+ if(!(last1pos == null)) last2pos = last1pos;
+ if(!(last0pos == null)) last1pos = last0pos;
+ if(!(last == null)) last0pos = last;
+
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+
+ lastDet = hit.getSubdetector().getName();
+ lastlayer5x5x2 = id.getNeighbourIDs(1, 2, 2);
+
+ }
+
+
+
+ if (exam > 2) {
+ if(debug) System.out.println("Debug: showering is already started " + exam + " layers ago" );
+ break;
+ }
+ }
+
+ if(debug){
+ System.out.println("Debug: last4pos= " + last4pos.magnitude() + " (" +last4pos.x() + " " + last4pos.y() + " " + last4pos.z()+ ")" );
+ System.out.println("Debug: last3pos= " + last3pos.magnitude() + " (" +last3pos.x() + " " + last3pos.y() + " " + last3pos.z()+ ")");
+ System.out.println("Debug: last2pos= " + last2pos.magnitude() + " (" + last2pos.x() + " " + last2pos.y() + " " + last2pos.z()+ ") ");
+ System.out.println("Debug: last1pos= " + last1pos.magnitude() + " (" + last1pos.x() + " " + last1pos.y() + " " + last1pos.z()+ ")");
+ System.out.println("Debug: last0pos= " + last0pos.magnitude() + " (" + last0pos.x() + " " + last0pos.y() + " " + last0pos.z()+ ")");
+ }
+
+ List<Hep3Vector> lastTwoPositions = new ArrayList<Hep3Vector>();
+ if(last4pos != null ){
+ lastTwoPositions.add(last3pos);
+ lastTwoPositions.add(last4pos);
+ }
+ else if(last3pos != null){
+ lastTwoPositions.add(last2pos);
+ lastTwoPositions.add(last3pos);
+ }
+ else if(last2pos != null){
+ lastTwoPositions.add(last1pos);
+ lastTwoPositions.add(last2pos);
+ }
+ else {
+ lastTwoPositions.add(last0pos);
+ lastTwoPositions.add(last1pos);
+ }
+
+
+
+ return lastTwoPositions;
+ }
+
+
+
+ protected Hep3Vector lastPosition(CalorimeterHit hit, int s, Set<Long> allHitIDs){
+
+ int count3x3=0;
+ int count5x5=0;
+ int count7x7=0;
+ int count9x9=0;
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int lay = id.getLayer();
+ double x = id.getPosition()[0];
+ double y = id.getPosition()[1];
+ double z = id.getPosition()[2];
+ Hep3Vector v = new BasicHep3Vector(x,y,z);
+ int p = (int) (v.magnitude());
+
+ long[] neighbours3x3 = id.getNeighbourIDs(0, 1, 1);
+ long[] neighbours5x5 = id.getNeighbourIDs(0, 2, 2);
+ long[] neighbours7x7 = id.getNeighbourIDs(0, 3, 3);
+ long[] neighbours9x9 = id.getNeighbourIDs(0, 4, 4);
+ for (long neighbourID3x3 : neighbours3x3)
+ if (allHitIDs.contains(neighbourID3x3)) count3x3++;
+ for (long neighbourID5x5 : neighbours5x5)
+ if (allHitIDs.contains(neighbourID5x5)) count5x5++;
+ for (long neighbourID7x7 : neighbours7x7)
+ if (allHitIDs.contains(neighbourID7x7)) count7x7++;
+ for (long neighbourID9x9 : neighbours9x9)
+ if (allHitIDs.contains(neighbourID9x9)) count9x9++;
+
+
+
+
+ if ( count9x9 >=2 ) exam++ ;
+ else exam = 0;
+
+ if(debug){
+ System.out.println("Debug: "+s+" hit: pos=" + p + " , " + count3x3 + " hits in 3x3, " + count5x5 + " hits in 5x5 " + count7x7 + " hits in 7x7 " + count9x9 + " hits in 9x9 " + lay + " layer of " + hit.getSubdetector().getName());
+ }
+
+ return v;
+ }
+
+
+}
