lcsim/src/org/lcsim/contrib/uiowa/MuonFinder
diff -N MipTrackMapDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MipTrackMapDriver.java 2 Dec 2008 00:15:42 -0000 1.1
@@ -0,0 +1,571 @@
+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.*;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * 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 MipTrackMapDriver extends Driver {
+ protected AIDA aida = AIDA.defaultInstance();
+ 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;
+
+ protected boolean use_DTree = false;
+
+ public MipTrackMapDriver (HelixExtrapolator extrap) {
+ m_extrap = extrap;
+ initTrackMatch();
+ }
+
+ public Map<Track, BasicCluster> createCalMIPMap(List<ReconstructedParticle> recoList) {
+ Map<Track, BasicCluster> MapTrkToMIP = new HashMap<Track, BasicCluster>();
+ for(ReconstructedParticle reco : recoList){
+ if(reco.getCharge() == 0) continue;
+ if(reco.getTracks().size() != 1) continue;
+ Track tr = reco.getTracks().get(0);
+ Set<CalorimeterHit> allhits = new HashSet<CalorimeterHit>();
+ for (Cluster recoClus : reco.getClusters()) allhits.addAll(recoClus.getCalorimeterHits());
+ m_allhits = allhits;
+ List<Cluster> allMatchableClusters = createClusters(allhits);
+ 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, 3, 3);
+ for (long lastnearID : nearbyHitArray ) { nearIDs.add(lastnearID); }
+
+ Hep3Vector last0 = new BasicHep3Vector(hit.getPosition());
+ Hep3Vector last1 = new BasicHep3Vector();
+ if( mipclus.getCalorimeterHits().size() > 1){
+ CalorimeterHit hit1 = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-2);
+ last1 = new BasicHep3Vector(hit1.getPosition());
+ }
+
+ //Searh cur hit in cur layer and arrange in angular order
+ //System.out.println("last = " + last0.x() + " " + last0.y() + " " + last0.z());
+ 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 InSide = (curlayer > lastlayer) && cursubdetName.contains(lastsubdetName) && nearIDs.contains(curhit.getCellID());
+ boolean EcalToHcalBarrel= (curid.getLayer()==0 && cursubdetName.contains("HADBarrel")) && lastsubdetName.contains("EMBarrel"); // ECal Barrel -> HCal Barrel
+ boolean EcalToHcalEndcap= (curid.getLayer()==0 && cursubdetName.contains("HADEndcap")) && lastsubdetName.contains("EMEndcap"); // ECal Endcap -> HCal Endcap
+ boolean BarrelEcalToEndHcal = (curid.getLayer()==0 && cursubdetName.contains("HADEndcap")) && lastsubdetName.contains("EMBarrel"); // ECal Barrel -> HCal Endcap (This is no way to HCal Barrel directly from ECal Endcap
+ boolean EndToBarrelEcal = (curid.getLayer()==0 && cursubdetName.contains("EMBarrel")) && lastsubdetName.contains("EMEndcap"); // ECal Endcap -> ECal Barrel
+ boolean EndToBarrelHcal = (curid.getLayer()==0 && cursubdetName.contains("HADBarrel")) && lastsubdetName.contains("HADEndcap"); // HCal Endcap -> HCal Barrel
+ if(EcalToHcalBarrel || EcalToHcalEndcap || InSide || EndToBarrelEcal || EndToBarrelHcal || BarrelEcalToEndHcal){
+ Hep3Vector curFromTangent = result.getTangent(curpos);
+ Hep3Vector curFromLastTwo = VecOp.sub(curpos, last0);
+ Hep3Vector cur = null;
+ if(curFromTangent != null) {
+ cur = curFromTangent;
+ }else {
+ cur = curFromLastTwo;
+ System.out.println("No getTangent From current");
+ }
+ Hep3Vector curUnit = VecOp.unit(cur);
+
+ Hep3Vector lastFromTangent = result.getTangent(last0);
+ Hep3Vector lastFromLastTwo = null;
+ if (mipclus.getCalorimeterHits().size() > 1) {
+ lastFromLastTwo = 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) {
+ lastFromLastTwo = result.getTangent();
+ }
+ if (lastFromLastTwo == null) {
+ // No extrapolation info and only one hit -- fall back to using
+ // a straight line from the origin.
+ lastFromLastTwo = VecOp.sub(last0, new BasicHep3Vector(0,0,0));
+
+ }
+ }
+ Hep3Vector last = null;
+ if(lastFromTangent != null) {
+ last = lastFromTangent;
+ }else {
+ last = lastFromLastTwo;
+ System.out.println("No getTangent From last");
+ }
+
+ Hep3Vector lastUnit = VecOp.unit(last);
+
+ //System.out.println("cur = " + curpos.x() + " " + curpos.y() + " " + curpos.z());
+ double dot = VecOp.dot(lastUnit, curUnit);
+ if(dot > 1) dot = 1.0;
+ double a = Math.acos(dot);
+ double d = VecOp.sub(curpos,last0).magnitude();
+ //System.out.println("a= " + a + " d= " +d);
+ 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 > 4) {
+ if(debug) System.out.println("Debug: showering is already started " + exam + " layers ago" );
+ double numRemove = 0;
+ if(mipclus.getCalorimeterHits().size() < 6) {
+ numRemove = mipclus.getCalorimeterHits().size() - 1 ;
+ } else numRemove = 5;
+ 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() < 2) k = mipclus.getCalorimeterHits().size();
+ else k = 2;
+ for(int i =0; i < k ;i++){
+ CalorimeterHit hit = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-1-i);
+ Hep3Vector v = new BasicHep3Vector(hit.getPosition());
+ System.out.println("Debug: last "+i+" pos= "+v.magnitude()+" (" +v.x() + " " + v.y() + " " + v.z()+ ")" );
+ }
+ System.out.println("cal MIP size= " + mipclus.getCalorimeterHits().size());
+ }
+ MapTrkToMIP.put(tr, mipclus);
+ } //track loop
+
+ return MapTrkToMIP;
+ }
+
+ public List<Cluster> createMIPMuDet(Set<CalorimeterHit> muhits){
+ removeUpperSensitiveLayer(muhits);
+ hitsInTree = new HitMap(muhits);
+ List<Cluster> muonmip = new Vector<Cluster>();
+
+ if(!use_DTree){
+ Set<Long> allHitIDs = new HashSet<Long>();
+ Set<CalorimeterHit> allhits = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> removedhits = new HashSet<CalorimeterHit>();
+ for(CalorimeterHit hit : muhits){
+ if(removedhits.contains(hit)) continue;
+ SortedMap<Double, Set<CalorimeterHit>> sortedHitbyPos= new TreeMap();
+ IDDecoder id = hit.getIDDecoder();
+ double layer = id.getLayer();
+ 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;
+ IDDecoder subid = subhit.getIDDecoder();
+ double sublayer = subid.getLayer();
+ Hep3Vector subhitPos = new BasicHep3Vector(subhit.getPosition());
+ Hep3Vector subhitPosUnit = VecOp.unit(subhitPos);
+ double cos = VecOp.dot(hitPosUnit, subhitPosUnit);
+ sortedHitbyAng.put(cos,subhit);
+ if(debug){
+ double cosfromIP = Math.abs(hitPos.z()/hitPos.magnitude());
+ if(cosfromIP < 0.8) {
+ aida.cloud1D("muon/cos between barrel hits from IP").fill(cos);
+ }else {
+ aida.cloud1D("muon/cos between endcap hits from IP").fill(cos);
+ }
+ }
+ 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);
+ }
+ }
+ muonmip.add(tmpClus);
+ }
+ }else{
+ List<Cluster> output = new Vector<Cluster>();
+ // 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);
+
+ System.out.println("size of output= " + output.size());
+ // Simple way by looking at direction of two clusters
+ List<Cluster> removed = new Vector<Cluster>();
+ if(output.size() > 1){
+ for(int i = 0 ; i < output.size() ; i ++){
+ Cluster cluster1 = output.get(i);
+ if(removed.contains(cluster1)) continue;
+ Hep3Vector v1 = getClusterDirection(output.get(i));
+ BasicCluster tmpClus = new BasicCluster();
+ for(CalorimeterHit hit : cluster1.getCalorimeterHits()){
+ tmpClus.addHit(hit);
+ }
+
+ for(int j = 1 ; j < output.size()-1 ; j ++){
+ Cluster cluster2 = output.get(j);
+ if(removed.contains(cluster2)) continue;
+ Hep3Vector v2 = getClusterDirection(output.get(j));
+ double cos = -1;
+ if( v1 != null && v2 !=null){ cos = VecOp.dot(v1, v2); }
+ if(cos>0.8) {
+ for(CalorimeterHit hit : cluster2.getCalorimeterHits()){
+ tmpClus.addHit(hit);
+ }
+ removed.add(cluster2);
+ }
+ }
+ muonmip.add(tmpClus);
+ removed.add(cluster1);
+ }
+ }else { muonmip.addAll(output); }
+ }
+
+ if(debug){
+ int hitsize = 0 ;
+ for(Cluster clus : muonmip){
+ hitsize = hitsize + clus.getCalorimeterHits().size();
+ for(CalorimeterHit hit : clus.getCalorimeterHits()){
+ examPosition(hit);
+ }
+ }
+ }
+
+ return muonmip;
+ }
+
+
+ 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);
+ System.out.println("Removed due to poor quality");
+ }
+ }
+ 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());
+ String detname = hit.getSubdetector().getName();
+
+ 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){
+ if(detname.contains("MuonEndcap") || detname.contains("MuonBarrel")){
+ long[] updownneighbours = id.getNeighbourIDs(1,0,0);
+ long[] sideneighbours = id.getNeighbourIDs(0,1,1);
+ for(long neighbourID : updownneighbours){
+ id.setID(neighbourID);
+ Hep3Vector v = id.getPositionVector();
+ double r = Math.sqrt(v.x()*v.x() + v.y()*v.y());
+ double z = v.z();
+ System.out.println("updown= " + neighbourID + " r= " + r + " z= " + z);
+ }
+ for(long neighbourID : sideneighbours){
+ id.setID(neighbourID);
+ Hep3Vector v = id.getPositionVector();
+ double r = Math.sqrt(v.x()*v.x() + v.y()*v.y());
+ double z = v.z();
+ System.out.println("side= " + neighbourID + " r= " + r + " z= " + z);
+ }
+ }
+ }
+
+ if(debug){
+ if(detname.contains("MuonEndcap") || detname.contains("MuonBarrel")){
+ Hep3Vector cellv = hitCellPosition(hit);
+ Hep3Vector v = new BasicHep3Vector(hit.getPosition());
+ int p = (int) (v.magnitude());
+ int r = (int) Math.sqrt(v.x()*v.x() + v.y()*v.y());
+ int z = (int) v.z();
+ int cellp = (int) (cellv.magnitude());
+ int cellr = (int) Math.sqrt(cellv.x()*cellv.x() + cellv.y()*cellv.y());
+ int cellz = (int) cellv.z();
+ System.out.println(" hit : p= " + p +" r=" + r + " , z " + z + ": cellID= " + hit.getCellID() + " " + count[0] + " hits in 3x3, " + count[1] + " at " + id.getLayer() + " layer of " + hit.getSubdetector().getName());
+ System.out.println(" cell: p= " + cellp +" r=" + cellr + " , z " + cellz);
+ }
+ }
+
+ }
+
+ protected void removeUpperSensitiveLayer(Set<CalorimeterHit> muhits){
+ Set<CalorimeterHit> upperhits = new HashSet<CalorimeterHit>();
+ for(CalorimeterHit hit : muhits){
+ String subdetName = hit.getSubdetector().getName();
+ Hep3Vector cellv = hitCellPosition(hit);
+ Hep3Vector v = new BasicHep3Vector(hit.getPosition());
+ int p = (int) (v.magnitude());
+ int r = (int) Math.sqrt(v.x()*v.x() + v.y()*v.y());
+ int z = (int) v.z();
+ int cellp = (int) (cellv.magnitude());
+ int cellr = (int) Math.sqrt(cellv.x()*cellv.x() + cellv.y()*cellv.y());
+ int cellz = (int) cellv.z();
+ boolean upper = true;
+ if(subdetName.contains("MuonEndcap")) upper = Math.abs(z) > Math.abs(cellz);
+ if(subdetName.contains("MuonBarrel")) upper = r > cellr;
+ if(upper) upperhits.add(hit);
+ }
+ muhits.removeAll(upperhits);
+ }
+
+ //give the vector of hit
+ 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;
+ }
+
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+
+ protected Hep3Vector getClusterDirection(Cluster clus){
+ Hep3Vector trackMuonUnit = new BasicHep3Vector();
+ if( clus.getCalorimeterHits().size() > 1){
+ CalorimeterHit muonhit1 = clus.getCalorimeterHits().get(clus.getCalorimeterHits().size()-1);
+ CalorimeterHit muonhit0 = clus.getCalorimeterHits().get(0);
+ Hep3Vector muonpos1 = new BasicHep3Vector(muonhit1.getPosition());
+ Hep3Vector muonpos0 = new BasicHep3Vector(muonhit0.getPosition());
+ Hep3Vector trackMuon = new BasicHep3Vector();
+ if(muonpos0.magnitude() > muonpos1.magnitude()){
+ trackMuon = VecOp.sub(muonpos0, muonpos1);
+ }else{
+ trackMuon = VecOp.sub(muonpos1, muonpos0);
+ }
+ trackMuonUnit = VecOp.unit(trackMuon);
+ } else { trackMuonUnit = null;}
+
+ return trackMuonUnit;
+ }
+
+ protected Set<Cluster> recursivelyFindSubClusters(Cluster clus) {
+ Set<Cluster> output = new HashSet<Cluster>();
+ output.add(clus);
+ for (Cluster subClus : clus.getClusters()) {
+ Set<Cluster> subClusterDaughters = recursivelyFindSubClusters(subClus);
+ output.addAll(subClusterDaughters);
+ }
+ return output;
+ }
+
+}
lcsim/src/org/lcsim/contrib/uiowa/MuonFinder
diff -u -r1.3 -r1.4
--- MuonFinderDriver.java 26 Nov 2008 23:21:27 -0000 1.3
+++ MuonFinderDriver.java 2 Dec 2008 00:15:42 -0000 1.4
@@ -22,9 +22,15 @@
import org.lcsim.geometry.Detector;
import org.lcsim.util.lcio.LCIOConstants;
import org.lcsim.util.aida.AIDA;
+import hep.physics.particle.Particle;
+import hep.physics.particle.properties.*;
+import org.lcsim.util.swim.HelixSwimmer;
+
/**
* Driver to find muon cluster and corresponding track.
- *
+ *
+ * This driver is only for reconstructed sample.
+ * This driver uses the reconstructed output.
* 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,
@@ -33,75 +39,65 @@
public class MuonFinderDriver extends Driver{
protected AIDA aida = AIDA.defaultInstance();
- protected String _tracklist;
- protected String _muonlist;
+ protected String _muonmap;
+ protected String _recolist;
protected double _bestmatch = 0.8;
- protected boolean _debug = false;
- protected boolean _debug1 = false;
+ protected boolean _debug = false;
protected boolean _useFirstTwoLayer = true;
- EventHeader m_event;
+ EventHeader _event;
//protected HelixExtrapolator _findCluster = new org.lcsim.recon.pfa.identifier.TrackHelixPlusHitExtrapolator();
protected HelixExtrapolator _findCluster = new org.lcsim.recon.pfa.identifier.TrackHelixExtrapolator();
+ //This extrapolator is only for standalone moun driver.
+
+ public MuonFinderDriver(String muonmap, String recolist){
+ _muonmap = muonmap;
+ _recolist = recolist;
+ // Filter muon system hits
+ {
+ DecisionMakerSingle<CalorimeterHit> upperLayer = new UpperSubLayerDecision();
+ DecisionMakerSingle<CalorimeterHit> lowerLayer = new NotDecisionMakerSingle(upperLayer);
+ DecisionMakerSingle<CalorimeterHit> timeCut = new CalorimeterHitTimeCutDecision(100);
+ AndDecisionMakerSingle<CalorimeterHit> lowerLayerAndTimeCut = new AndDecisionMakerSingle<CalorimeterHit>();
+ lowerLayerAndTimeCut.addDecisionMaker(lowerLayer);
+ lowerLayerAndTimeCut.addDecisionMaker(timeCut);
+ add(new ListFilterDriver(lowerLayerAndTimeCut, "MuonBarrHits", "CorrMuonBarrDigiHits", CalorimeterHit.class));
+ add(new ListFilterDriver(lowerLayerAndTimeCut, "MuonEndcapHits", "CorrMuonEndcapDigiHits", CalorimeterHit.class));
+ }
- public MuonFinderDriver(String tracklist, String muonlist){
- _tracklist = tracklist;
- _muonlist = muonlist;
}
public void process(EventHeader event)
{
super.process(event);
- m_event = event;
+ _event = event;
- Set<CalorimeterHit> calhits = new HashSet<CalorimeterHit>();
+ List<Track> trackList = new Vector<Track>();
+ List<ReconstructedParticle> recoParticles=event.get(ReconstructedParticle.class, _recolist);
+ for(ReconstructedParticle reco : recoParticles){
+ ParticleType type = ParticlePropertyManager.getParticlePropertyProvider().get(13);
+ BaseParticleID pid = new BaseParticleID(type);
+ double currentParticleMass = type.getMass();
+
+ if(reco.getCharge() == 0) continue;
+ if(reco.getTracks().size() != 1) continue;
+ Track recoTrack = reco.getTracks().get(0);
+ trackList.add(recoTrack);
+ }
Set<CalorimeterHit> mudethits = new HashSet<CalorimeterHit>();
-
- List<CalorimeterHit> allHitsEcalBarrel = event.get(CalorimeterHit.class, "EcalBarrDigiHits");
- List<CalorimeterHit> allHitsEcalEndcap = event.get(CalorimeterHit.class, "EcalEndcapDigiHits");
- List<CalorimeterHit> allHitsHcalBarrel = event.get(CalorimeterHit.class, "HcalBarrDigiHits");
- List<CalorimeterHit> allHitsHcalEndcap = event.get(CalorimeterHit.class, "HcalEndcapDigiHits");
- calhits.addAll(allHitsEcalBarrel);
- calhits.addAll(allHitsEcalEndcap);
- calhits.addAll(allHitsHcalBarrel);
- calhits.addAll(allHitsHcalEndcap);
List<CalorimeterHit> allHitsMUBarrel = event.get(CalorimeterHit.class, "CorrMuonBarrDigiHits");
List<CalorimeterHit> allHitsMUEndcap = event.get(CalorimeterHit.class, "CorrMuonEndcapDigiHits");
mudethits.addAll(allHitsMUBarrel);
mudethits.addAll(allHitsMUEndcap);
- HitMap calHitMap = new HitMap();
- HitMap muHitMap = new HitMap();
- for (CalorimeterHit hit : calhits) {
- long id = hit.getCellID();
- calHitMap.put(id, hit);
- }
-
- for (CalorimeterHit hit : mudethits) {
- long id = hit.getCellID();
- muHitMap.put(id, hit);
- }
-
- Set<Long> allHitIDs = new HashSet<Long>();
- for (CalorimeterHit allhit : calhits) allHitIDs.add(allhit.getCellID());
-
_findCluster.process(event); // picks up geometry
- List<Track> trackList = event.get(Track.class, _tracklist);
-
- //Create mip cluster separately in Calorimetry and Muon detector.
- //Mip in CAL has track matched
- //But Mip in Muon detector is standalone.
- //In the future, the muon Mip will be replace with standalone muon track.
- MipTrackMap finder = new MipTrackMap(_findCluster);
- Map<Track,BasicCluster> mipmap = finder.createCalMIPMap(calhits, trackList);
+ MipTrackMapDriver finder = new MipTrackMapDriver(_findCluster);
+ Map<Track,BasicCluster> mipmap = finder.createCalMIPMap(recoParticles);
List<Cluster> muonmips = finder.createMIPMuDet(mudethits);
Map<Track,Set<Cluster>> outputmap = new HashMap<Track,Set<Cluster>>();
- if(_debug){
- System.out.println("Muon mip size= " + muonmips.size()+" (should be 1 for signle sample)");
- aida.cloud1D("muon/muon mip size").fill(muonmips.size());
- }
+
for(Cluster mumip : muonmips){
if(_debug) System.out.println("Candidate Muon");
if( mumip.getCalorimeterHits().size() > 1){
@@ -168,7 +164,7 @@
Hep3Vector hitpos = new BasicHep3Vector();
Hep3Vector lastpoint = new BasicHep3Vector();
- boolean muon = false;
+ //boolean muon = false;
int isolated = 0;
if(calmip == null || result == null) {
if(_debug) System.out.println(" REASON: Null mip or extrapolation");
@@ -180,7 +176,8 @@
int layer = id.getLayer();
String subdetName = hit.getSubdetector().getName();
//Followings are required to take pion track out for muon candidate
- muon = subdetName.contains("HADBarrel") || (subdetName.contains("HADEndcap") && (layer > 29) );
+ //muon = subdetName.contains("HADBarrel") || subdetName.contains("HADEndcap");
+ //This requirement is not used for this standalone muon driver
hitpos = new BasicHep3Vector(hit.getPosition());
//Using extrapolated track
int lastLayerBarrel = getLastLayer("HADBarrel");
@@ -204,10 +201,8 @@
tpoint = tEndcap;
isBarrel = false;
}else if(tBarrel != null && tEndcap != null){
- //System.out.println("Both have it");
tpoint = tBarrel;
}else{
- //System.out.println("Both don't have it");
for(int l=0; l<endingLayer+1 ; l++){
int nextlayer = endingLayer-l;
if(isBarrel) {tpoint = result.extendToHCALBarrelLayer(nextlayer);}
@@ -226,25 +221,16 @@
else {tid = event.getDetector().getDecoder("HcalEndcapHits");}
long cellID = tid.findCellContainingXYZ(tpoint);
tid.setID(cellID);
- long[] neighbours = tid.getNeighbourIDs(0,5,5);
- int count = 0;
- if(calHitMap.containsKey(cellID)) {count++;}
- for(long neighbour : neighbours){
- if(calHitMap.containsKey(neighbour)){count++;}
- }
-
if(i==0){
lastpoint = tpoint;
last = result.getTangent(tpoint);
lastUnit = VecOp.unit(last);
+ break; // Only for Stand alone muon driver to find out JUST the last point
}
- double rpos = Math.sqrt(tpoint.x()*tpoint.x() + tpoint.y()*tpoint.y());
- if(_debug1) System.out.println("Extrapolated track at " + endingLayer + " : r= "+rpos+" z= "+tpoint.z() + " hit= " + count );
if(endingLayer == 0){
endingLayer = 40;
isBarrel = false;
}
- if(count > 0 && count < 3) isolated++;
}
if(!hasLastPoint){
if(_debug) System.out.println(" REASON: Null extrapolated track last point");
@@ -255,18 +241,10 @@
continue;
}
- //Require last position or 7 isolated layer in last 10 layers obtained by extrapolated track.
- //This helps remove pion which is showering in CAL and also take the MIP which is contaminated by
- //other shower so it stopped in the middle of CAL.
- //if(isolated > 6) {
- // if(_debug) System.out.println("Isolation cut failed");
- // continue;
- //}
- if (!muon) {
- // Skip these cases to avoid pion contamination.
- if(_debug) System.out.println(" REASON: Mip doesn't reach in last 10 layers");
- continue;
- }
+ //The requirement for the last mip point used to be here but
+ //remove the last mip point requirement
+ //since this can get better efficiency for stand alone muon driver
+
Hep3Vector link = VecOp.sub(muonpos0, lastpoint);
Hep3Vector linkUnit = VecOp.unit(link);
double cos0 = VecOp.dot(lastUnit, MuonUnit);
@@ -334,9 +312,11 @@
List<Track> muonTracks = new Vector<Track>(MuonMapTrackToCluster.keySet());
List<Cluster> muonClusters = new Vector<Cluster>(MuonMapTrackToCluster.values());
int flag = 1<<LCIOConstants.CLBIT_HITS;
- event.put(_muonlist, muonTracks, Track.class, flag);
- if(_debug) event.put("MuonClusters", muonClusters, Cluster.class, flag );
- if(_debug) System.out.println("Muon size= " + muonTracks.size());
+ event.put(_muonmap, MuonMapTrackToCluster);
+ if(_debug){
+ event.put("Muons", muonTracks, Track.class, flag);
+ event.put("MuonClusters", muonClusters, Cluster.class, flag );
+ }
}
public Hep3Vector getHitMeanPosition(Collection<CalorimeterHit> hits){
@@ -352,9 +332,9 @@
}
public int getLastLayer(String s){
- Detector det = m_event.getDetector();
+ Detector det = _event.getDetector();
CylindricalCalorimeter detname = ((CylindricalCalorimeter) det.getSubdetectors().get(s));
return detname.getLayering().getLayerCount();
}
-
+
}