lcsim/src/org/lcsim/contrib/uiowa/MuonFinder
diff -N MipTrackMap.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MipTrackMap.java 28 Sep 2008 06:31:15 -0000 1.1
@@ -0,0 +1,365 @@
+package org.lcsim.contrib.uiowa.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.*;
+
+/**
+ * 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.
+ */
+
+
+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;
+
+
+
+ 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, 2, 2);
+ for (long lastnearID : nearbyHitArray ) { nearIDs.add(lastnearID); }
+
+ Hep3Vector last0 = hitPosition(hit);
+ Hep3Vector last1 = new BasicHep3Vector();
+ if( mipclus.getCalorimeterHits().size() > 1){
+ CalorimeterHit hit1 = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-2);
+ last1 = hitPosition(hit1);
+ }
+
+ //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 EcalToHcal= lastlayer == 30 && curp > p && curid.getLayer() == 0;
+ boolean HcalToMuDet = lastlayer == 33 && curp > p && curid.getLayer() == 0;
+ boolean InSide = (curp > p || curlayer > lastlayer) && nearIDs.contains(curhit.getCellID());
+ boolean EndToBarrelEcal = (curid.getLayer()==0 && cursubdetName.contains("EMBarrel")) && lastsubdetName.contains("EMEndcap");
+ boolean EndToBarrelHcal = (curid.getLayer()==0 && cursubdetName.contains("HADBarrel")) && lastsubdetName.contains("HADEndcap");
+ if(debug){
+ if(EndToBarrelEcal || EndToBarrelHcal) System.out.println("This is a boundary condition");
+ }
+
+ if(EcalToHcal || InSide || EndToBarrelEcal || EndToBarrelHcal ){
+ Hep3Vector cur = VecOp.sub(curpos, last0);
+ Hep3Vector curUnit = VecOp.unit(cur);
+ Hep3Vector last = null;
+ if (mipclus.getCalorimeterHits().size() > 1) {
+ last = 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) {
+ last = result.getTangent();
+ }
+ if (last == null) {
+ // No extrapolation info and only one hit -- fall back to using
+ // a straight line from the origin.
+ last = VecOp.sub(last0, new BasicHep3Vector(0,0,0));
+ }
+ }
+ Hep3Vector lastUnit = VecOp.unit(last);
+ double a = Math.acos(VecOp.dot(lastUnit, curUnit));
+ double d = VecOp.sub(cur,last).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 > 2) {
+ if(debug) System.out.println("Debug: showering is already started " + exam + " layers ago" );
+ double numRemove = 0;
+ if(mipclus.getCalorimeterHits().size() < 4) {
+ numRemove = mipclus.getCalorimeterHits().size() - 1 ;
+ } else numRemove = 3;
+ 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() < 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;
+ }
+
+ public List<Cluster> createMIPMuDet(Set<CalorimeterHit> muhits){
+ List<Cluster> output = new Vector<Cluster>();
+
+ hitsInTree = new HitMap(muhits);
+ // 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 (m_removePoorQualityMips) {
+ removePoorQualityMips(output);
+ }
+
+ return output;
+ }
+
+
+ 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);
+ }
+ }
+ 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 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;
+ }
+
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+
+}