lcsim/src/org/lcsim/recon/cluster/mipfinder
diff -N InitialMipFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ InitialMipFinder.java 6 Oct 2010 20:01:26 -0000 1.1
@@ -0,0 +1,237 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.cluster.mipfinder;
+import org.lcsim.recon.cluster.mipfinder.trackxtrap.*;
+import org.lcsim.event.CalorimeterHit;
+import java.util.*;
+import org.lcsim.recon.cluster.util.*;
+import hep.physics.vec.*;
+import org.lcsim.geometry.Calorimeter;
+import org.lcsim.util.hitmap.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class InitialMipFinder
+{
+ protected TrackXtrapCalPositions txp;
+ CalorimeterHitDSort chds;
+ int[] multhitinlaycut = {2,2};
+ int[] nconsecmhitcut = {3,3};
+ int[] nlaymisscut = {3,3};
+ double[] dsqcutv;
+ double[] dsqcuta;
+ double[] ncelldcutv = {1.5,2.5};
+ double[] ncelldcuta = {2.5,3.5};
+ //TODO - find the cellsize
+ double[] cellsize = {3.5,10.};
+ CalorimeterHit[] bestHits;
+ CalorimeterHit[][] allCutHits;
+ int[] nallPerStep;
+ int[] ncutPerStep;
+ int[] indPerStep;
+ boolean onePerLayer = false;
+ int maxhitsperlayerinmip = 5;
+ Hep3Vector[] off2;
+ public InitialMipFinder()
+ {
+ chds = new CalorimeterHitDSort();
+ dsqcutv = new double[cellsize.length];
+ dsqcuta = new double[cellsize.length];
+ for(int i=0;i<cellsize.length;i++)
+ {
+ dsqcutv[i] = cellsize[i]*cellsize[i]*ncelldcutv[i]*ncelldcutv[i];
+ dsqcuta[i] = cellsize[i]*cellsize[i]*ncelldcuta[i]*ncelldcuta[i];
+ }
+ }
+ public void setOnePerLayer(boolean x){onePerLayer = x;}
+ public void setMultHitInLayCut(int[] x){multhitinlaycut = x;}
+ public void setNConsecMHitCut(int[] x){nconsecmhitcut = x;}
+ public void setNLayMissCut(int[] x){nlaymisscut = x;}
+ public void setNCellDCutV(double[] x)
+ {
+ ncelldcutv = x;
+ for(int i=0;i<cellsize.length;i++)
+ {
+ dsqcutv[i] = cellsize[i]*cellsize[i]*ncelldcutv[i]*ncelldcutv[i];
+ }
+ }
+ public void setNCellDCutA(double[] x)
+ {
+ ncelldcuta = x;
+ for(int i=0;i<cellsize.length;i++)
+ {
+ dsqcuta[i] = cellsize[i]*cellsize[i]*ncelldcuta[i]*ncelldcuta[i];
+ }
+ }
+ public BasicCluster findMip(TrackXtrapInfo txi, HitMap hitmap)
+ {
+ bestHits = findMip0(txi,hitmap);
+ if(bestHits == null)
+ {
+ System.out.println("No hits in early layer within cut");
+ return null;
+ }
+ BasicCluster mipcl = new BasicCluster();
+ int ncon0 = 0;
+ int nconm = 0;
+ int lasti = bestHits.length - 1;
+// String vstr = "Valid hits: ";
+// String astr = "Nearby hits: ";
+// for(int i=0;i<bestHits.length;i++)
+// {
+// vstr += ncutPerStep[i]+", ";
+// astr += nallPerStep[i]+", ";
+// }
+// System.out.println(vstr);
+// System.out.println(astr);
+ for(int i=0;i<bestHits.length;i++)
+ {
+ if(ncutPerStep[i] == 0)ncon0++;
+ else ncon0 = 0;
+ if(ncon0 >= nlaymisscut[indPerStep[i]])
+ {
+ lasti = i - nlaymisscut[indPerStep[i]];
+ break;
+ }
+ if(nallPerStep[i] >= multhitinlaycut[indPerStep[i]])nconm++;
+ else nconm = 0;
+ if(nconm >= nconsecmhitcut[indPerStep[i]])
+ {
+ lasti = i - nconsecmhitcut[indPerStep[i]];
+ if(lasti < 0)lasti = 0;
+ break;
+ }
+ }
+// System.out.println("Last valid step = "+lasti);
+ if(lasti < 0)return null;
+ if(onePerLayer)
+ {
+ for(int i=0;i<lasti+1;i++)
+ {
+ if(bestHits[i] != null)mipcl.addHit(bestHits[i]);
+ }
+ }
+ else
+ {
+ for(int i=0;i<lasti+1;i++)
+ {
+ if(ncutPerStep[i] > maxhitsperlayerinmip)
+ {
+ if(bestHits[i] != null)mipcl.addHit(bestHits[i]);
+ }
+ else
+ {
+ for(int j=0;j<Math.min(20,ncutPerStep[i]);j++)
+ {
+ mipcl.addHit(allCutHits[i][j]);
+ }
+ }
+ }
+ }
+ if(mipcl.getCalorimeterHits().size() > 0)return mipcl;
+ return null;
+ }
+ public CalorimeterHit[] findMip0(TrackXtrapInfo txi, HitMap hmap)
+ {
+ List<XtrapPosition> xpl = txi.getPositions();
+ if(xpl.size() < 1)return null;
+ CalorimeterHit[] bhits = new CalorimeterHit[xpl.size()];
+ Hep3Vector[] offsets = new Hep3Vector[xpl.size()];
+ off2 = new Hep3Vector[xpl.size()];
+ nallPerStep = new int[xpl.size()];
+ ncutPerStep = new int[xpl.size()];
+ allCutHits = new CalorimeterHit[xpl.size()][20];
+ indPerStep = new int[xpl.size()];
+ for(int i=0;i<xpl.size();i++)
+ {
+ XtrapPosition xp = xpl.get(i);
+ int ic = 0;
+ if( (xp.getCalorimeterType() == Calorimeter.CalorimeterType.HAD_BARREL)||
+ (xp.getCalorimeterType() == Calorimeter.CalorimeterType.HAD_ENDCAP) )ic = 1;
+ indPerStep[i] = ic;
+ List<CalorimeterHit> hl = txi.getHits(xp);
+ Hep3Vector p = xp.getPosition();
+ if(hl.size() < 1)
+ {
+ nallPerStep[i] = 0;
+ ncutPerStep[i] = 0;
+ bhits[i] = null;
+ continue;
+ }
+ Hep3Vector ap = p;
+ if(i > 20)ap = adjustPosition(i,p,offsets);
+ CalorimeterHit bestHit = null;
+ double bestdsq = 9999.;
+ int ncut = 0;
+ int nall = 0;
+ for(CalorimeterHit h:hl)
+ {
+ if(hmap.containsValue(h))
+ {
+ double dsq = getDsqToHit(ap,h);
+ if(dsq < dsqcutv[ic])
+ {
+ if(ncut < 20)allCutHits[i][ncut] = h;
+ ncut++;
+ if(dsq < bestdsq)
+ {
+ bestdsq = dsq;
+ bestHit = h;
+ }
+ }
+ if(dsq < dsqcuta[ic])
+ {
+ nall++;
+ }
+ }
+ }
+ bhits[i] = bestHit;
+ if(bestHit != null)
+ {
+ offsets[i] = VecOp.sub(new BasicHep3Vector(bestHit.getPosition()), p);
+ off2[i] = VecOp.sub(new BasicHep3Vector(bestHit.getPosition()), ap);
+ }
+ ncutPerStep[i] = ncut;
+ nallPerStep[i] = nall;
+ }
+ return bhits;
+ }
+ protected Hep3Vector adjustPosition(int i,Hep3Vector p,Hep3Vector[] off)
+ {
+ int j = i-1;
+ int n = 0;
+ Hep3Vector totoff = new BasicHep3Vector();
+ while( (j >= 0)&&(n < 5))
+ {
+ if(off[j] != null)
+ {
+ totoff = VecOp.add(totoff, off[j]);
+ n++;
+ }
+ j--;
+ }
+ if(n != 5)return p;
+ totoff = VecOp.mult(.2, totoff);
+ return VecOp.add(p,totoff);
+ }
+ protected double getDsqToHit(Hep3Vector pos,CalorimeterHit h)
+ {
+ Hep3Vector ihp = new BasicHep3Vector(h.getPosition());
+ return VecOp.sub(pos, ihp).magnitudeSquared();
+ }
+ public int[] getNAllPerStep()
+ {
+ return nallPerStep;
+ }
+ public int[] getNValidPerStep()
+ {
+ return ncutPerStep;
+ }
+ public CalorimeterHit[][] getAllValidHits(){return allCutHits;}
+ public Hep3Vector[] getOffsets(){return off2;}
+}
lcsim/src/org/lcsim/recon/cluster/mipfinder
diff -N ShowerPointFinderDriver2.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ShowerPointFinderDriver2.java 6 Oct 2010 20:01:27 -0000 1.1
@@ -0,0 +1,94 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.cluster.mipfinder;
+import org.lcsim.recon.cluster.mipfinder.trackxtrap.*;
+import java.util.*;
+import org.lcsim.util.*;
+import org.lcsim.event.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.recon.cluster.util.*;
+
+/**
+ * Wrapper class for ShowerPointFinder2, allowing it to be
+ * run as a Driver. Outputs are written to the event.
+ * Assumes track extrapolations and potential hit assignments have already been done.
+ *
+ * @author cassell
+ */
+public class ShowerPointFinderDriver2 extends Driver
+{
+ protected String m_inputTrackCalhitMap;
+ protected String m_inputTrackXtrapInfo;
+ protected String m_inputHitMap;
+ protected String m_inputTrackList;
+ protected String m_outputTrackMap;
+ protected String m_outputHitMap;
+ protected String m_outputMipClusterList;
+ protected boolean first;
+ InitialMipFinder imp;
+ /**
+ * Constructor. Quite a few inputs needed:
+ *
+ * @param inputHitMap Name of the input hitmap.
+ * @param inputTrackXtrapInfo Name of the input TrackXtrapInfo list.
+ * @param outputTrackMap Name to write an output Map<Track,BasicCluster> out as, connecting tracks to their MIP stubs.
+ * @param outputHitMap Name to write the remaining, unused hits out as.
+ * @param outputMipClusterList Name to write the MIP stub clusters out as.
+ */
+ public ShowerPointFinderDriver2(String inputHitMap, String inputTrackList, String inputTrackXtrapInfo, String outputTrackMap, String outputHitMap, String outputMipClusterList)
+ {
+ m_inputHitMap = inputHitMap;
+ m_inputTrackList = inputTrackList;
+ m_outputTrackMap = outputTrackMap;
+ m_outputHitMap = outputHitMap;
+ m_outputMipClusterList = outputMipClusterList;
+ m_inputTrackXtrapInfo = inputTrackXtrapInfo;
+ first = true;
+ }
+
+ public void process(EventHeader event) {
+ if(first)
+ {
+ first = false;
+ imp = new InitialMipFinder();
+ }
+ super.process(event);
+ // Read in & set up
+ HitMap inputHitMap = (HitMap)(event.get(m_inputHitMap));
+ HitMap outputHitMap = new HitMap(inputHitMap);
+ List<Track> inputTracks = event.get(Track.class, m_inputTrackList);
+ Map<Track,BasicCluster> mapTrkToMIP = new HashMap<Track,BasicCluster>();
+ // Loop over tracks and find initial mip.
+ List<TrackXtrapInfo> txil = new ArrayList<TrackXtrapInfo>(event.get(TrackXtrapInfo.class,m_inputTrackXtrapInfo));
+ for(TrackXtrapInfo txi:txil)
+ {
+ Track t = txi.getTrack();
+ if(inputTracks.contains(t))
+ {
+ BasicCluster c = imp.findMip(txi,inputHitMap);
+ if(c != null)mapTrkToMIP.put(t,c);
+ }
+ }
+ // Remaining, unused hits
+ for (BasicCluster clus : mapTrkToMIP.values()) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ long id = hit.getCellID();
+ outputHitMap.remove(id);
+ }
+ }
+
+ // Write main output
+ event.put(m_outputTrackMap, mapTrkToMIP);
+ event.put(m_outputHitMap, outputHitMap);
+
+ // Additional list of MIP clusters used to identify MIP hits for dE/dx etc:
+ List<Cluster> preShowerMips = new Vector<Cluster>();
+ preShowerMips.addAll(mapTrkToMIP.values());
+ event.put(m_outputMipClusterList, preShowerMips);
+ event.getMetaData(preShowerMips).setTransient(true);
+ }
+
+}
