lcsim/src/org/lcsim/recon/cluster/muonfinder
diff -N MuonFinderWrapper3.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonFinderWrapper3.java 6 Oct 2010 20:00:22 -0000 1.1
@@ -0,0 +1,118 @@
+package org.lcsim.recon.cluster.muonfinder;
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.util.*;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.HitInMCALDecision;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.recon.util.CalorimeterInformation;
+import org.lcsim.geometry.Calorimeter.CalorimeterType;
+import org.lcsim.recon.cluster.mipfinder.trackxtrap.TrackXtrapInfo;
+
+/**
+ * Wrapper class for MuonFinder2.
+ *
+ * @author [log in to unmask]
+ */
+
+public class MuonFinderWrapper3 extends Driver {
+
+ // Inputs
+ protected String m_trackList;
+ protected String m_inHitMap;
+ protected String m_trackXtrapInfo;
+
+ // Outputs
+ protected String m_outMuonTrackClusterMap;
+ protected String m_outHitMap;
+ protected String m_outTrackList;
+ protected String _inMu = "tmpMu";
+ protected MuonFinder3 mf2;
+ protected int nSensitivePerLayer;
+ protected boolean init;
+ CalorimeterInformation ci;
+
+ public MuonFinderWrapper3(String inTrackXtrapInfo,String inTrackList, String inHitMap, String outMuonTrackClusterMap, String outHitMap, String outTrackList)
+ {
+ // Inputs/outputs
+// System.out.println("Initializing MuonFinderWrapper2");
+ m_trackXtrapInfo = inTrackXtrapInfo;
+ m_trackList = inTrackList;
+ m_inHitMap = inHitMap;
+ m_outMuonTrackClusterMap = outMuonTrackClusterMap;
+ m_outHitMap = outHitMap;
+ m_outTrackList = outTrackList;
+ mf2 = new MuonFinder3();
+ init = false;
+ }
+
+ public void process(EventHeader event) {
+ if(!init)
+ {
+ init = true;
+ ci = CalorimeterInformation.instance();
+ nSensitivePerLayer = ci.getSubdetector(CalorimeterType.MUON_BARREL).getLayering().getLayer(0).getSensorIndices().size();
+ if(event.getDetectorName().contains("sid02"))nSensitivePerLayer = 1;
+ }
+ // Read in
+ HitMap inHitMap = (HitMap)(event.get(m_inHitMap));
+ List<CalorimeterHit> inHitList = new ArrayList<CalorimeterHit>(inHitMap.values());
+ List<Track> inTrackList = event.get(Track.class, m_trackList);
+
+ // Massage inputs into the form wanted by old process routine
+ HitMap muonHits = new HitMap();
+ HitMap otherHits = new HitMap();
+ HitInMCALDecision dec = new HitInMCALDecision();
+ for (CalorimeterHit hit : inHitMap.values()) {
+ long id = hit.getCellID();
+ if (dec.valid(hit)) {
+ muonHits.put(id, hit);
+ }
+ else otherHits.put(id, hit);
+ }
+
+ // Prepare outputs
+ List<Track> outTrackList = new Vector<Track>();
+ HitMap outHitMap = new HitMap(inHitMap);
+
+ // Do the work
+ Map<Track,Cluster> outMuonTrackClusterMap = mf2.findMuons(nSensitivePerLayer, otherHits, new ArrayList<CalorimeterHit>(muonHits.values()),inTrackList,event.get(TrackXtrapInfo.class,m_trackXtrapInfo));
+
+ // Make HitMap of remaining hits:
+ for (Cluster muonClus : outMuonTrackClusterMap.values()) {
+ for (CalorimeterHit hit : muonClus.getCalorimeterHits()) {
+ long id = hit.getCellID();
+ outHitMap.remove(id);
+ }
+ }
+
+ // Make List of remaining tracks:
+ outTrackList.addAll(inTrackList);
+ outTrackList.removeAll(outMuonTrackClusterMap.keySet());
+
+
+ // Write out
+ event.put(m_outHitMap, outHitMap);
+ event.put(m_outTrackList, outTrackList);
+
+ // One last thing: We occasionally get two muon clusters that overlap.
+ // Go through and prevent that:
+ Map<Track,Cluster> nonOverlappingOutputMap = new HashMap<Track,Cluster>();
+ Set<Long> usedHits = new HashSet<Long>();
+ for (Track muonTrack : outMuonTrackClusterMap.keySet()) {
+ Cluster oldMuonCluster = outMuonTrackClusterMap.get(muonTrack);
+ BasicCluster newMuonCluster = new BasicCluster();
+ for (CalorimeterHit hit : oldMuonCluster.getCalorimeterHits()) {
+ long id = hit.getCellID();
+ if (!usedHits.contains(id)) {
+ usedHits.add(id);
+ newMuonCluster.addHit(hit);
+ }
+ nonOverlappingOutputMap.put(muonTrack, newMuonCluster);
+ }
+ }
+ event.put(m_outMuonTrackClusterMap, nonOverlappingOutputMap);
+ }
+}
+
lcsim/src/org/lcsim/recon/cluster/muonfinder
diff -N MuonFinder3.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonFinder3.java 6 Oct 2010 20:00:22 -0000 1.1
@@ -0,0 +1,313 @@
+package org.lcsim.recon.cluster.muonfinder;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.*;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.Track;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.MuonCalorimeterHitLSort;
+import org.lcsim.recon.util.CalorimeterInformation;
+import org.lcsim.geometry.Calorimeter.CalorimeterType;
+import org.lcsim.geometry.Calorimeter;
+import org.lcsim.recon.cluster.mipfinder.trackxtrap.TrackXtrapInfo;
+import org.lcsim.recon.cluster.mipfinder.InitialMipFinder;
+import org.lcsim.util.hitmap.HitMap;
+
+
+/**
+ * Driver to find muon cluster and corresponding track.
+ *
+ * search mip cluster in Muon detector and assume that that cluster is from muon
+ * then find track matched with this mip cluster by looking at direction.
+ * Since pion can reach muon detector as well specially in endcap,
+ * it is required to find good quality mip cluster in muon detector.
+ *
+ * @author [log in to unmask]
+ * @version $Id: MuonFinder3.java,v 1.1 2010/10/06 20:00:22 cassell Exp $
+ */
+
+public class MuonFinder3 {
+ protected double ctcut = 0.85;
+ protected boolean _useFirstTwoLayer = true;
+
+ protected boolean init;
+ CalorimeterInformation ci;
+ Comparator<CalorimeterHit> chc;
+ double tpcut = 3.;
+ InitialMipFinder imp;
+ int[] mhilc = {3,3};
+ int[] ncmhc = {5,5};
+ int[] mslc = {4,4};
+ double[] ncdcv = {3.5,6.};
+ double[] ncdca = {5.,6.};
+ int maxconnoniso = 20;
+
+
+ public MuonFinder3(){
+ init = false;
+ chc = new MuonCalorimeterHitLSort();
+ imp = new InitialMipFinder();
+ imp.setMultHitInLayCut(mhilc);
+ imp.setNConsecMHitCut(mhilc);
+ imp.setNLayMissCut(mslc);
+ imp.setNCellDCutV(ncdcv);
+ imp.setNCellDCutA(ncdca);
+ }
+
+ public Map<Track,Cluster> findMuons(int nspl, HitMap otherHits, List<CalorimeterHit> muhits, List<Track> tracklist, List<TrackXtrapInfo> txilist)
+ {
+ Map<Track,Cluster> outMap = new HashMap<Track,Cluster>();
+ if(!init)
+ {
+ init = true;
+ ci = CalorimeterInformation.instance();
+ }
+
+ // Find mips in the muon detectors
+ List<Cluster> mumips = findMuonMips(muhits,nspl);
+ if(mumips.size() < 1)
+ {
+ return outMap;
+ }
+ Map<Cluster,TrackXtrapInfo> matches = new HashMap<Cluster,TrackXtrapInfo>();
+ // Loop over mumips
+ for(Cluster mumip:mumips)
+ {
+ // Sort the hits by distance from origin
+ List<CalorimeterHit> mhl = mumip.getCalorimeterHits();
+ Collections.sort(mhl,chc);
+ // Use hits in the first 2 layers to get a direction
+ List<CalorimeterHit> flhits = new ArrayList<CalorimeterHit>();
+ List<CalorimeterHit> slhits = new ArrayList<CalorimeterHit>();
+ CalorimeterHit fhit = mhl.get(0);
+ IDDecoder fid = fhit.getIDDecoder();
+ fid.setID(fhit.getCellID());
+ int flay = fid.getVLayer();
+ int slay = -1;
+ CalorimeterType sct = null;
+ CalorimeterType fct = ((Calorimeter)(fid.getSubdetector())).getCalorimeterType();
+ flhits.add(fhit);
+ int nextind = 1;
+ for(int i=1;i<mhl.size();i++)
+ {
+ CalorimeterHit hit = mhl.get(i);
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int lay = id.getVLayer();
+ CalorimeterType ct = ((Calorimeter)(id.getSubdetector())).getCalorimeterType();
+ if( (lay == flay)&&(ct == fct) )
+ {
+ flhits.add(hit);
+ }
+ else
+ {
+ slay = lay;
+ sct = ct;
+ nextind = i;
+ slhits.add(hit);
+ break;
+ }
+ }
+ for(int i=nextind+1;i<mhl.size();i++)
+ {
+ CalorimeterHit hit = mhl.get(i);
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int lay = id.getVLayer();
+ CalorimeterType ct = ((Calorimeter)(id.getSubdetector())).getCalorimeterType();
+ if( (lay == slay)&&(ct == sct) )slhits.add(hit);
+ else break;
+ }
+ double pos0x = 0.;
+ double pos0y = 0.;
+ double pos0z = 0.;
+ for(CalorimeterHit h:flhits)
+ {
+ double[] pos = h.getPosition();
+ pos0x += pos[0];
+ pos0y += pos[1];
+ pos0z += pos[2];
+ }
+ int nh = flhits.size();
+ Hep3Vector fpos = new BasicHep3Vector(pos0x/nh, pos0y/nh, pos0z/nh);
+ double pos1x = 0.;
+ double pos1y = 0.;
+ double pos1z = 0.;
+ for(CalorimeterHit h:slhits)
+ {
+ double[] pos = h.getPosition();
+ pos1x += pos[0];
+ pos1y += pos[1];
+ pos1z += pos[2];
+ }
+ nh = slhits.size();
+ Hep3Vector spos = new BasicHep3Vector(pos1x/nh, pos1y/nh, pos1z/nh);
+ Hep3Vector mudir = VecOp.sub(spos, fpos);
+ Hep3Vector muunit = VecOp.unit(mudir);
+ // Match mumips with xtrapolated tracks
+ TrackXtrapInfo bestmatch = null;
+ double bestavct = -1.;
+ int nmm = 0;
+ for(TrackXtrapInfo txi:txilist)
+ {
+ Track t = txi.getTrack();
+ if(tracklist.contains(t))
+ {
+ if((new BasicHep3Vector(t.getMomentum())).magnitude() > tpcut)
+ {
+ Hep3Vector tpos = txi.getPositions().get(txi.getPositions().size() -1).getPosition();
+ Hep3Vector tunit = txi.getPositions().get(txi.getPositions().size() -1).getDirection();
+ Hep3Vector tmudir = VecOp.sub(fpos,tpos);
+ Hep3Vector tmuunit = VecOp.unit(tmudir);
+ // As per TaeJeong, average 3 combinations of 3 directions.
+ double ct0 = VecOp.dot(tunit, tmuunit);
+ double ct1 = VecOp.dot(tunit, muunit);
+ double ct2 = VecOp.dot(tmuunit, muunit);
+ double cta = (ct0+ct1+ct2)/3.;
+ if(cta > bestavct)
+ {
+ bestmatch = txi;
+ bestavct = cta;
+ }
+ }
+ }
+ }
+ if(bestmatch != null)
+ {
+ if(bestavct > ctcut)
+ {
+ matches.put(mumip, bestmatch);
+ }
+ else
+ {
+ nmm++;
+ continue;
+ }
+ }
+ }
+ // Now do the mip finding in the calorimeters
+ for(Cluster mucl:matches.keySet())
+ {
+ TrackXtrapInfo txi = matches.get(mucl);
+ CalorimeterHit[] hits = imp.findMip0(txi,otherHits);
+ // Cut out some poorly matched calorimeter mips
+ // As per TaeJeong, require > 6 isolated hits in last 20 layers
+ int[] nallperstep = imp.getNAllPerStep();
+ int[] nvperstep = imp.getNValidPerStep();
+ // If less than 50 crossings with a hit, reject
+ if (nallperstep.length < 50)continue;
+ int niso = 0;
+ for(int i = nallperstep.length-1;i>nallperstep.length-20;i--)
+ {
+ if( (nvperstep[i] > 0)&&(nallperstep[i] < 3) )niso++;
+ }
+ if(niso < 7)
+ {
+ continue;
+ }
+ // Add addition requirement for no more than maxconnoniso
+ // consecutive nonisolated hits.
+ int nisoa = 0;
+ int ncnon = 0;
+ int maxncnon = 0;
+ for(int i=0;i<nallperstep.length;i++)
+ {
+ if( (nvperstep[i] > 0)&&(nallperstep[i] < 3) )
+ {
+ nisoa++;
+ if(ncnon > maxncnon)maxncnon = ncnon;
+ ncnon = 0;
+ }
+ else ncnon++;
+ }
+ if(ncnon > maxncnon)maxncnon = ncnon;
+ if(maxncnon > maxconnoniso)continue;
+ // Make the cluster and call this a muon
+ BasicCluster mcl = new BasicCluster();
+ mcl.addCluster(mucl);
+ for(int i=0;i<hits.length;i++)
+ {
+ if(hits[i] != null)mcl.addHit(hits[i]);
+ }
+ // Make the output map
+ outMap.put(txi.getTrack(), mcl);
+ }
+ return outMap;
+ }
+ public List<Cluster> findMuonMips(List<CalorimeterHit> muhits,int napl)
+ {
+ List<Cluster> outList = new ArrayList<Cluster>();
+ if(napl == 1)removeUpperSensitiveLayer(muhits);
+ Set<CalorimeterHit> removedhits = new HashSet<CalorimeterHit>();
+ for(CalorimeterHit hit : muhits){
+ if(removedhits.contains(hit)) continue;
+ SortedMap<Double, Set<CalorimeterHit>> sortedHitbyPos= new TreeMap();
+ Hep3Vector hitPos = new BasicHep3Vector(hit.getPosition());
+ Hep3Vector hitPosUnit = VecOp.unit(hitPos);
+
+ Set<CalorimeterHit> hits = sortedHitbyPos.get(hitPos.magnitude());
+ if(hits == null){
+ hits = new HashSet<CalorimeterHit>();
+ sortedHitbyPos.put(hitPos.magnitude(),hits);
+ }
+ hits.add(hit);
+
+ removedhits.add(hit);
+ SortedMap<Double, CalorimeterHit> sortedHitbyAng= new TreeMap();
+ for(CalorimeterHit subhit : muhits){
+ if(removedhits.contains(subhit)) continue;
+ Hep3Vector subhitPos = new BasicHep3Vector(subhit.getPosition());
+ Hep3Vector subhitPosUnit = VecOp.unit(subhitPos);
+ double cos = VecOp.dot(hitPosUnit, subhitPosUnit);
+ sortedHitbyAng.put(cos,subhit);
+ if(cos > 0.96) {
+
+ Set<CalorimeterHit> subhits = sortedHitbyPos.get(subhitPos.magnitude());
+ if(subhits == null){
+ subhits = new HashSet<CalorimeterHit>();
+ sortedHitbyPos.put(subhitPos.magnitude(),subhits);
+ }
+ if(subhit == null) { System.out.println("ERROR no subhit"); }
+ subhits.add(subhit);
+ removedhits.add(subhit);
+ }
+ }
+ BasicCluster tmpClus = new BasicCluster();
+ for(Set<CalorimeterHit> c : sortedHitbyPos.values()){
+ for(CalorimeterHit h : c){
+ tmpClus.addHit(h);
+ }
+ }
+ SimpleMipQualityDecision check = new SimpleMipQualityDecision(napl);
+ if(check.valid(tmpClus))outList.add(tmpClus);
+ }
+ return outList;
+ }
+ protected void removeUpperSensitiveLayer(List<CalorimeterHit> muhits){
+ Set<CalorimeterHit> upperhits = new HashSet<CalorimeterHit>();
+ double err = .01;
+ for(CalorimeterHit hit : muhits){
+ Hep3Vector cellv = hitCellPosition(hit);
+ Hep3Vector v = new BasicHep3Vector(hit.getPosition());
+ double delta = v.magnitude() - cellv.magnitude();
+ if(delta > err)upperhits.add(hit);
+ }
+ muhits.removeAll(upperhits);
+ }
+ protected Hep3Vector hitCellPosition(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;
+ }
+
+}