lcsim/src/org/lcsim/recon/cluster/muonfinder
diff -N MuonDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonDriver.java 10 Dec 2008 19:54:09 -0000 1.1
@@ -0,0 +1,315 @@
+package org.lcsim.recon.cluster.muonfinder;
+
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.util.*;
+import org.lcsim.event.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.base.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.recon.cluster.directedtree.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+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.geometry.IDDecoder;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+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,
+ * it is required to find good quality mip cluster in muon detector.
+ */
+
+public class MuonDriver extends Driver{
+ protected AIDA aida = AIDA.defaultInstance();
+ protected String _muonmap;
+ protected String _recolist;
+ protected double _bestmatch = 0.8;
+ protected boolean _debug = false;
+ protected boolean _useFirstTwoLayer = true;
+ 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 MuonDriver(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 void process(EventHeader event)
+ {
+ super.process(event);
+ _event = event;
+
+ List<Track> trackList = new Vector<Track>();
+ List<ReconstructedParticle> recoParticles=event.get(ReconstructedParticle.class, _recolist);
+ for(ReconstructedParticle reco : recoParticles){
+ ParticleType type = ParticlePropertyManager.getParticlePropertyProvider().get(11);
+ BaseParticleID pid = new BaseParticleID(type);
+ double electronMass = type.getMass();
+
+ if(reco.getCharge() == 0) continue;
+ if(reco.getMass() == electronMass) continue;
+ if(reco.getTracks().size() != 1) continue;
+ Track recoTrack = reco.getTracks().get(0);
+ trackList.add(recoTrack);
+ }
+ Set<CalorimeterHit> mudethits = new HashSet<CalorimeterHit>();
+ List<CalorimeterHit> allHitsMUBarrel = event.get(CalorimeterHit.class, "CorrMuonBarrDigiHits");
+ List<CalorimeterHit> allHitsMUEndcap = event.get(CalorimeterHit.class, "CorrMuonEndcapDigiHits");
+ mudethits.addAll(allHitsMUBarrel);
+ mudethits.addAll(allHitsMUEndcap);
+
+ _findCluster.process(event); // picks up geometry
+
+ MipTrackMap finder = new MipTrackMap(_findCluster);
+ List<Cluster> muonmips = finder.createMIPMuDet(mudethits);
+
+ Map<Track,Set<Cluster>> outputmap = new HashMap<Track,Set<Cluster>>();
+
+ for(Cluster mumip : muonmips){
+ if(_debug) System.out.println("Candidate Muon");
+ if( mumip.getCalorimeterHits().size() > 1){
+ Hep3Vector muonpos0 = new BasicHep3Vector();
+ Hep3Vector muonpos1 = new BasicHep3Vector();
+ //Muon detector has 20cm iron between 2 layers. The direction is not affected by the segmentation.
+ //Trying to use first two layers to identify direction intead of using first and last hit.
+ //In endcap, track curves so it will bring wrong direction if we use first and last hit.
+ //First two layers will give better direction.
+ if(_useFirstTwoLayer){
+ SortedMap<Integer, Set<CalorimeterHit>> exam = new TreeMap();
+ int prelayer = -1 ;
+ int i = -1;
+ String predetName = null;
+ for(CalorimeterHit h : mumip.getCalorimeterHits()){
+ IDDecoder id = h.getIDDecoder();
+ id.setID(h.getCellID());
+ int layer = id.getLayer();
+ if(prelayer != layer || prelayer == -1){
+ predetName = h.getSubdetector().getName();
+ i++;
+ if(i==2) break;
+ }
+
+ Set<CalorimeterHit> hits = exam.get(layer);
+ if(hits == null){
+ hits = new HashSet<CalorimeterHit>();
+ exam.put(layer,hits);
+ }
+ String subdetName = h.getSubdetector().getName();
+ if(predetName.contains(subdetName)) hits.add(h);
+ prelayer = layer;
+ }
+
+ Collection<CalorimeterHit> hits0 = exam.get(exam.firstKey());
+ Collection<CalorimeterHit> hits1 = exam.get(exam.lastKey());
+ muonpos0 = getHitMeanPosition(hits0);
+ muonpos1 = getHitMeanPosition(hits1);
+
+ }else{
+ CalorimeterHit hit0 = mumip.getCalorimeterHits().get(0);
+ CalorimeterHit hit1 = mumip.getCalorimeterHits().get(mumip.getCalorimeterHits().size()-1);
+ muonpos0 = new BasicHep3Vector(hit0.getPosition());
+ muonpos1 = new BasicHep3Vector(hit1.getPosition());
+ }
+
+ Hep3Vector Muon = VecOp.sub(muonpos1, muonpos0);
+ Hep3Vector MuonUnit = VecOp.unit(Muon);
+
+ double bestmatch = -1;
+ Track bestTrack = null;
+ Cluster bestmipc = null;
+ for(Track tr : trackList){
+ if(_debug) System.out.println(" Candidate Track");
+ Hep3Vector track = (new BasicHep3Vector(tr.getMomentum()));
+ double p = track.magnitude();
+ double pT = Math.sqrt(track.x()*track.x() + track.y()*track.y());
+ double costheta = Math.abs(track.z()/p);
+
+ HelixExtrapolationResult result = _findCluster.performExtrapolation(tr);
+ Hep3Vector last = new BasicHep3Vector();
+ Hep3Vector lastUnit = new BasicHep3Vector();
+ Hep3Vector lastpoint = new BasicHep3Vector();
+
+ if(p < 3) continue; //It would not work below 3 GeV
+ if(result == null) {
+ if(_debug) System.out.println(" DEBUG: Null mip or extrapolation");
+ continue;
+ }else {
+ //Using extrapolated track
+ int lastLayerBarrel = getLastLayer("HADBarrel");
+ int lastLayerEndcap = getLastLayer("HADEndcap");
+ if(lastLayerBarrel != lastLayerEndcap) {System.out.println("ERROR : number of layer different");}
+ int endingLayer = lastLayerBarrel;
+ boolean isBarrel = true;
+ boolean hasLastPoint = true;
+ for(int i=0 ; i < 20 ;i++){
+ endingLayer = endingLayer - 1;
+ IDDecoder tid = null;
+ Hep3Vector tpoint = new BasicHep3Vector();
+ Hep3Vector tBarrel = null;
+ if(isBarrel) {
+ tBarrel = result.extendToHCALBarrelLayer(endingLayer);
+ }
+ Hep3Vector tEndcap = result.extendToHCALEndcapLayer(endingLayer);
+ if(tBarrel != null && tEndcap == null){
+ tpoint = tBarrel;
+ }else if(tBarrel == null && tEndcap != null){
+ tpoint = tEndcap;
+ isBarrel = false;
+ }else if(tBarrel != null && tEndcap != null){
+ tpoint = tBarrel;
+ }else{
+ for(int l=0; l<endingLayer+1 ; l++){
+ int nextlayer = endingLayer-l;
+ if(isBarrel) {tpoint = result.extendToHCALBarrelLayer(nextlayer);}
+ else {tpoint = result.extendToHCALEndcapLayer(nextlayer);}
+ if(tpoint != null){
+ endingLayer = nextlayer;
+ break;
+ }
+ }
+ }
+ if(tpoint == null) {
+ hasLastPoint = false;
+ break;
+ }
+ if(isBarrel){tid = event.getDetector().getDecoder("HcalBarrHits");}
+ else {tid = event.getDetector().getDecoder("HcalEndcapHits");}
+ long cellID = tid.findCellContainingXYZ(tpoint);
+ tid.setID(cellID);
+ 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
+ }
+ if(endingLayer == 0){
+ endingLayer = 40;
+ isBarrel = false;
+ }
+ }
+ if(!hasLastPoint){
+ if(_debug) System.out.println(" DEBUG: Null extrapolated track last point");
+ 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);
+ double cos1 = VecOp.dot(lastUnit, linkUnit);
+ double cos2 = VecOp.dot(MuonUnit, linkUnit);
+ double cos = (cos0 + cos1 + cos2)/3;
+ //find best matched track
+ //compare the combination of three directions
+ if( cos > bestmatch ) {
+ if(_debug){
+ aida.cloud1D("muon/cos_extrapolatedtrack").fill(cos0);
+ aida.cloud1D("muon/cos_extrapolatedtrack_link").fill(cos1);
+ aida.cloud1D("muon/cos_muon_link").fill(cos2);
+ aida.cloud1D("muon/combined cos").fill(cos);
+ System.out.println(" cos0= " + cos0 + " cos1= " + cos1 + " cos2= " + cos2 + " average cos= " + cos);
+ System.out.println(" DEBUG: This is possible MUON Track!");
+ }
+ bestmatch = cos;
+ bestTrack = tr;
+ }else { if(_debug) System.out.println(" DEBUG: This is NOT possible muon Track!"); }
+ }
+ if(bestmatch > _bestmatch){
+ Set<Cluster> c = outputmap.get(bestTrack);
+ if (c == null) {
+ c = new HashSet<Cluster>();
+ outputmap.put(bestTrack, c);
+ }
+ c.add(mumip);
+ if(_debug) System.out.println("DEBUG: Found one track for this muon mip");
+ }else {
+ if(_debug) {
+ System.out.println("DEBUG: This cluster in MuDet has no track matched!" + " due to " + bestmatch);
+ }
+ }
+ }else { if(_debug) System.out.println ("DEBUG: This muon mip has 0 number of hits!");}
+ }
+
+ Map<Track, Cluster> MuonMapTrackToCluster = new HashMap<Track,Cluster>();
+ for(Track tr : outputmap.keySet()){
+ BasicCluster tmpClus = new BasicCluster();
+ for(Cluster clus : outputmap.get(tr)){
+ for(CalorimeterHit hit : clus.getCalorimeterHits()){
+ tmpClus.addHit(hit);
+ }
+ }
+
+ //Check Quality for mip in Muon Detector only
+ SimpleMipQualityDecision check = new SimpleMipQualityDecision();
+ boolean passesCheck = check.valid(tmpClus);
+ if(passesCheck) {
+ MuonMapTrackToCluster.put(tr,tmpClus);
+ }
+ }
+
+ List<Track> muonTracks = new Vector<Track>(MuonMapTrackToCluster.keySet());
+ List<Cluster> muonClusters = new Vector<Cluster>(MuonMapTrackToCluster.values());
+ int flag = 1<<LCIOConstants.CLBIT_HITS;
+ 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){
+ Hep3Vector hitPosition = new BasicHep3Vector();
+ double size = hits.size();
+ double norm = 1/size;
+ for(CalorimeterHit hit : hits) {
+ hitPosition = VecOp.add(hitPosition, new BasicHep3Vector(hit.getPosition()));
+ }
+ Hep3Vector meanPosition = VecOp.mult(norm, hitPosition);
+
+ return meanPosition;
+ }
+
+ public int getLastLayer(String s){
+ Detector det = _event.getDetector();
+ CylindricalCalorimeter detname = ((CylindricalCalorimeter) det.getSubdetectors().get(s));
+ return detname.getLayering().getLayerCount();
+ }
+
+}