Commit in lcsim/src/org/lcsim/recon/cluster/density on MAIN
DTreeAnalysis.java+170-1511.2 -> 1.3
DirectedTreeClusterer.java+59-381.2 -> 1.3
+229-189
2 modified files
GL: Add energy threshold, reduce output and cleanup

lcsim/src/org/lcsim/recon/cluster/density
DTreeAnalysis.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- DTreeAnalysis.java	12 Dec 2005 05:15:17 -0000	1.2
+++ DTreeAnalysis.java	27 Dec 2005 16:33:50 -0000	1.3
@@ -13,15 +13,22 @@
 import org.lcsim.event.CalorimeterHit;
 import org.lcsim.event.MCParticle;
 import org.lcsim.geometry.compact.Readout;
-import org.lcsim.geometry.segmentation.SegmentationBase;
 import org.lcsim.digisim.CellSelector;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
 import org.lcsim.recon.cluster.util.BasicCluster;
 import org.lcsim.recon.cluster.util.SortClustersBySize;
+import org.lcsim.geometry.IDDecoder;
 
 public class DTreeAnalysis extends Driver {
 
-    public DTreeAnalysis() throws IOException {
+    private int _debug=0;
+
+    public DTreeAnalysis() {
         _runPar = RunControlParameters.getInstance();
+	_emEcut = _runPar.getEMmip() * _runPar.getEMthresh() / _runPar.getEMweight();
+	_hdEcut = _runPar.getHDmip() * _runPar.getHDthresh() / _runPar.getHDweight();
+
+	// aida stuff
 	AIDA aida = AIDA.defaultInstance();
         _tree = aida.tree();
         _tf = aida.analysisFactory().createTupleFactory(_tree);
@@ -33,33 +40,37 @@
 
     public void process(EventHeader event) {
 
-	System.out.println("*************** New Event **********");
         _evtnum = event.getEventNumber();
-	System.out.println("Starting DTreeAnalysis for event no. "+_evtnum);
+	if(_debug>0) {
+	    System.out.println("Starting DTreeAnalysis, event # "+_evtnum);
+	}
         _tree.cd("DirectedTree");
 
 	String embName = "EcalBarrHits";
 	String emeName = "EcalEndcapHits";
 	String hdbName = "HcalBarrHits";
 	String hdeName = "HcalEndcapHits";
-	embhitmap = _loader.getCollHitMap(embName);
-	emehitmap = _loader.getCollHitMap(emeName);
-	hdbhitmap = _loader.getCollHitMap(hdbName);
-	hdehitmap = _loader.getCollHitMap(hdeName);
-
-	System.out.println("DTree: #hits: EMB="+embhitmap.size()
-			   +", EMEC="+emehitmap.size()
-			   +", HB="+hdbhitmap.size()
-			   +", HEC="+hdehitmap.size());
+	embhitmap = _expert.getCollHitMap(embName, _emEcut);
+	emehitmap = _expert.getCollHitMap(emeName, _emEcut);
+	hdbhitmap = _expert.getCollHitMap(hdbName, _hdEcut);
+	hdehitmap = _expert.getCollHitMap(hdeName, _hdEcut);
+
+	if(_debug>0) {
+	  System.out.println("DTree: #hits: "
+			     +" EMB="+embhitmap.size()
+			     +", EMEC="+emehitmap.size()
+			     +", HDB="+hdbhitmap.size()
+			     +", HDEC="+hdehitmap.size());
+	}
 
         if(_runPar.ClusterSeparately()){
             _calType = "EM";
-	    processCollection(event,embName);
-// 	    processCollection(event,emeName);
+ 	    processCollection(event,embName,embhitmap);
+//  	    processCollection(event,emeName,emehitmap);
 
             _calType = "HD";
-	    processCollection(event,hdbName);
-// 	    processCollection(event,hdeName);
+ 	    processCollection(event,hdbName,hdbhitmap);
+//  	    processCollection(event,hdeName,hdehitmap);
         }
         else {
 	  assert false : "Sorry, single-pass clustering unavailable for now.";
@@ -89,7 +100,7 @@
             Class[] columnClassesMCem = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
             _tupleMCem = _tf.create("MCEM","Gen Clusters in EM",columnNamesMCem,columnClassesMCem);
 
-            String[] columnNamesMChd = {"nclusMChd","hdmcClusFolder={float ene,float x,float y,float z,int size,int pid,int chrg,float e,float px,float py,float pz,float vx,float vy,float vz}","hdmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi,float r,float phi,float theta,float time}","evtno"};
+            String[] columnNamesMChd = {"nclusMChd","hdmcClusFolder={float ene,float x,float y,float z,int size,int pid,float chrg,float e,float px,float py,float pz,float vx,float vy,float vz}","hdmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi,float r,float phi,float theta,float time}","evtno"};
             Class[] columnClassesMChd = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
             _tupleMChd = _tf.create("MCHD","Gen Clusters in HD",columnNamesMChd,columnClassesMChd);
         }
@@ -98,7 +109,8 @@
     }
 
     public void fillDirectedTreeNtuple(List<BasicCluster> recon,
-				       List<BasicCluster> mctruth)
+				       List<BasicCluster> mctruth,
+				       IDDecoder decoder)
     {
         if(_calType=="EM"){
             emClusFolder = _tupleEM.getTuple(1);
@@ -114,48 +126,48 @@
 	    hdmcCellFolder = _tupleMChd.getTuple(2);
         }
 
+	// Fill tuples with reco info
         int tag = 0;
 	Map<CalorimeterHit,Integer> recTagMap
 	    = new HashMap<CalorimeterHit,Integer>();
 	for(BasicCluster clust : recon ) {
             List<CalorimeterHit> cellVec = clust.getCalorimeterHits();
-	    System.out.println(tag+" "+cellVec.size());
-	    if(cellVec.size()==0) System.out.println("*** bad cellVec in "+_calType);
-//             for(int i=0;i<cellVec.size();i++){
-//                 MyCalorimeterHit myhit = (MyCalorimeterHit) cellVec.get(i);
+// 	    System.out.println(tag+" "+cellVec.size());
+	    if(cellVec.size()==0)
+		System.out.println("*** bad cellVec in "+_calType);
 	    for(CalorimeterHit ihit : cellVec) {
 		Integer tagObject = new Integer(tag);
-// 		System.out.println("putting <"+myhit+"->"+tagObject+">");
+//  		System.out.println("putting <"+MyTools.printID(ihit.getCellID())+" -> "+tagObject+">");
 		recTagMap.put(ihit,tagObject);
                 double[] pos = ihit.getPosition();
 		double cdens = _loader.getDensity(ihit);
 		long jid = ihit.getCellID();
-		int ily = MyTools.getLayer(jid);
+		int ily1 = MyTools.getLayer(jid);
 		int iz = MyTools.getThetaBin(jid);
 		int iphi = MyTools.getPhiBin(jid);
-                if(_calType=="EM"){
-                    emCellFolder.fill(0,(float) ihit.getRawEnergy());
-                    emCellFolder.fill(1,(float)pos[0]);
-                    emCellFolder.fill(2,(float)pos[1]);
-                    emCellFolder.fill(3,(float)pos[2]);
-		    emCellFolder.fill(4,(float) cdens);
-                    emCellFolder.fill(5,tag);
-		    emCellFolder.fill(6,ily);
-		    emCellFolder.fill(7,iz);
-		    emCellFolder.fill(8,iphi);
-                    emCellFolder.addRow();
-                }
-                if(_calType=="HD"){
-                    hdCellFolder.fill(0,(float) ihit.getRawEnergy());
-                    hdCellFolder.fill(1,(float)pos[0]);
-                    hdCellFolder.fill(2,(float)pos[1]);
-                    hdCellFolder.fill(3,(float)pos[2]);
-		    hdCellFolder.fill(4,(float) cdens);
-                    hdCellFolder.fill(5,tag);
-                    hdCellFolder.fill(6,ily);
-                    hdCellFolder.fill(7,iz);
-                    hdCellFolder.fill(8,iphi);
-                    hdCellFolder.addRow();
+		decoder.setID( ihit.getCellID() );
+		int ily = decoder.getValue(0);
+		int iu = decoder.getValue(3);
+		int iv = decoder.getValue(4);
+// 		System.out.println("ilay="+ily1+" "+ily+", iu="+iphi+" "+iu+", iv="+iz+" "+iv);
+		assert ily==ily1 : "problem in layer index";
+		assert iphi==iu : "problem in phi/iU index";
+		assert iz==iv : "problem in iz/iV index";
+
+		ITuple cellFolder = null;
+		if(_calType=="EM") cellFolder = emCellFolder;
+		if(_calType=="HD") cellFolder = hdCellFolder;
+		if(cellFolder!=null) {
+                    cellFolder.fill(0,(float) ihit.getRawEnergy());
+                    cellFolder.fill(1,(float)pos[0]);
+                    cellFolder.fill(2,(float)pos[1]);
+                    cellFolder.fill(3,(float)pos[2]);
+		    cellFolder.fill(4,(float)cdens);
+                    cellFolder.fill(5,tag);
+		    cellFolder.fill(6,ily);
+		    cellFolder.fill(7,iu);
+		    cellFolder.fill(8,iv);
+                    cellFolder.addRow();
                 }
             }
 	    ITuple clusFolder = null;
@@ -163,24 +175,17 @@
 	    if(_calType=="HD") clusFolder = hdClusFolder;
 	    if(clusFolder!=null) {
 		Hep3Vector pos = new BasicHep3Vector( clust.getPosition() );
-                clusFolder.fill(0,(float) clust.getEnergy());
+                clusFolder.fill(0,(float) clust.getRawEnergy());
                 clusFolder.fill(1,(float) pos.x());
                 clusFolder.fill(2,(float) pos.y());
                 clusFolder.fill(3,(float) pos.z());
                 clusFolder.fill(4, clust.getSize());
                 clusFolder.addRow();
             }
-//             if(_calType=="HD"){
-//                 hdClusFolder.fill(0,(float) clust.getEnergy());
-//                 hdClusFolder.fill(1,(float) clust._xpos);
-//                 hdClusFolder.fill(2,(float) clust._ypos);
-//                 hdClusFolder.fill(3,(float) clust._zpos);
-//                 hdClusFolder.fill(4, clust.getSize() );
-//                 hdClusFolder.addRow();
-//             }
             tag++;
 	}
 
+	// Fill tuples with MC info
         tag = 0;
 	for( BasicCluster clust : mctruth ) {
 	  for( CalorimeterHit ihit : clust.getCalorimeterHits() ) {
@@ -194,9 +199,18 @@
 	    double cdens = _loader.getDensity(ihit);
 	    if(cdens==0) continue;  // skip hits below energy cut
 	    long jid = ihit.getCellID();
-	    int ily = MyTools.getLayer(jid);
+	    int ily1 = MyTools.getLayer(jid);
 	    int iz = MyTools.getThetaBin(jid);
 	    int iphi = MyTools.getPhiBin(jid);
+	    decoder.setID( ihit.getCellID() );
+	    int ily = decoder.getValue(0);
+	    int iu = decoder.getValue(3);
+	    int iv = decoder.getValue(4);
+// 	    System.out.println("ilay="+ily1+" "+ily+", iu="+iz+" "+iu+", iv="+iphi+" "+iv);
+	    assert ily==ily1 : "problem in layer index";
+	    assert iphi==iv : "problem in phi/iV index";
+	    assert iz==iu : "problem in iz/iU index";
+
 	    double cellR = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
 	    double theta = Math.atan2(cellR,pos[2]);
 	    double phi = Math.atan2(pos[1],pos[0]);
@@ -204,108 +218,82 @@
 	    SimCalorimeterHit iihit = (SimCalorimeterHit)ihit;
 	    int npart = iihit.getMCParticleCount();
 	    double timeStamp = 99999.;
-	    for(int j=0;j<npart;j++){
+	    for(int j=0;j<npart;j++) {
 		double tmpTime = (double) iihit.getContributedTime(j);
 		if(tmpTime<timeStamp){
 		    timeStamp = tmpTime;
 		}
 	    }
-	    if(_calType=="EM") {
-	      emmcCellFolder.fill(0,(float) ihit.getRawEnergy());
-	      emmcCellFolder.fill(1,(float)pos[0]);
-	      emmcCellFolder.fill(2,(float)pos[1]);
-	      emmcCellFolder.fill(3,(float)pos[2]);
-	      emmcCellFolder.fill(4,(float) cdens);
-	      emmcCellFolder.fill(5,tag);
-	      emmcCellFolder.fill(6,rtag);
-	      emmcCellFolder.fill(7,ily);
-	      emmcCellFolder.fill(8,iz);
-	      emmcCellFolder.fill(9,iphi);
-	      emmcCellFolder.fill(10,(float) cellR);
-	      emmcCellFolder.fill(11,(float) phi);
-	      emmcCellFolder.fill(12,(float) theta);
-	      emmcCellFolder.fill(13,(float) timeStamp);
-	      emmcCellFolder.addRow();
+
+	    ITuple mcCellFolder = null;
+	    if(_calType=="EM") mcCellFolder = emmcCellFolder;
+	    if(_calType=="HD") mcCellFolder = hdmcCellFolder;
+	    if(mcCellFolder!=null) {
+	      mcCellFolder.fill(0,(float) ihit.getRawEnergy());
+	      mcCellFolder.fill(1,(float)pos[0]);
+	      mcCellFolder.fill(2,(float)pos[1]);
+	      mcCellFolder.fill(3,(float)pos[2]);
+	      mcCellFolder.fill(4,(float) cdens);
+// 	      System.out.println("Gtag="+tag+", Rtag="+rtag);
+	      mcCellFolder.fill(5,tag);
+	      mcCellFolder.fill(6,rtag);
+	      mcCellFolder.fill(7,ily);
+	      mcCellFolder.fill(8,iu);
+	      mcCellFolder.fill(9,iv);
+	      mcCellFolder.fill(10,(float) cellR);
+	      mcCellFolder.fill(11,(float) phi);
+	      mcCellFolder.fill(12,(float) theta);
+	      mcCellFolder.fill(13,(float) timeStamp);
+	      mcCellFolder.addRow();
 //  	      System.out.println("EM Clus cells tag="+tag+", rtag="+rtag
 //  				 +", id="+MyTools.printID(jid)
 // // 				 +", pos="+pos[0]+" "+pos[1]+" "+pos[2]
 //  				 +", dens="+cdens);
 	    }
-	    if(_calType=="HD") {
-	      hdmcCellFolder.fill(0,(float) ihit.getRawEnergy());
-	      hdmcCellFolder.fill(1,(float)pos[0]);
-	      hdmcCellFolder.fill(2,(float)pos[1]);
-	      hdmcCellFolder.fill(3,(float)pos[2]);
-	      hdmcCellFolder.fill(4,(float) cdens);
-	      hdmcCellFolder.fill(5,tag);
-	      hdmcCellFolder.fill(6,rtag);
-	      hdmcCellFolder.fill(7,ily);
-	      hdmcCellFolder.fill(8,iz);
-	      hdmcCellFolder.fill(9,iphi);
-	      hdmcCellFolder.fill(10,(float) cellR);
-	      hdmcCellFolder.fill(11,(float) phi);
-	      hdmcCellFolder.fill(12,(float) theta);
-	      hdmcCellFolder.fill(13,(float) timeStamp);
-	      hdmcCellFolder.addRow();
-// 	      if(tag==7){
-// 		  System.out.println("HD Clus cells: id="+MyTools.printID(jid)
-// 				     +", pos="+pos[0]+" "+pos[1]+" "+pos[2]
-// 				     +", dens="+cdens);
-// 	      }
-	    }
 	  }
-	  MCParticle mcpart = getMCParticleInCluster(clust);
 
-	  float charge = (float)mcpart.getCharge();
-	  int particleID = mcpart.getType().getPDGID();
-	  float mcE = (float)mcpart.getEnergy();
-	  Hep3Vector mcMom = mcpart.getMomentum();
-	  double[] pos = clust.getPosition();
-	  Hep3Vector vert = mcpart.getOrigin();
-	  if(_calType=="EM"){
-	    emmcClusFolder.fill(0,(float) clust.getEnergy());
-	    emmcClusFolder.fill(1, (float)pos[0] );
-	    emmcClusFolder.fill(2, (float)pos[1] );
-	    emmcClusFolder.fill(3, (float)pos[2] );
-	    emmcClusFolder.fill(4, clust.getSize() );
-	    emmcClusFolder.fill(5, particleID);
-	    emmcClusFolder.fill(6, charge);
-	    emmcClusFolder.fill(7, mcE);
-	    emmcClusFolder.fill(8, (float)mcMom.x());
-	    emmcClusFolder.fill(9,(float)mcMom.y());
-	    emmcClusFolder.fill(10,(float)mcMom.z());
-	    emmcClusFolder.fill(11,(float) vert.x());
-	    emmcClusFolder.fill(12,(float) vert.y());
-	    emmcClusFolder.fill(13,(float) vert.z());
-	    emmcClusFolder.addRow();
-	  }
-	  if(_calType=="HD") {
-	    hdmcClusFolder.fill(0, (float)clust.getEnergy() );
-	    hdmcClusFolder.fill(1, (float)pos[0] );
-	    hdmcClusFolder.fill(2, (float)pos[1] );
-	    hdmcClusFolder.fill(3, (float)pos[2] );
-	    hdmcClusFolder.fill(4, clust.getSize() );
-	    hdmcClusFolder.fill(5, particleID);
-	    hdmcClusFolder.fill(6, (int) charge);
-	    hdmcClusFolder.fill(7, mcE);
-	    hdmcClusFolder.fill(8, (float)mcMom.x());
-	    hdmcClusFolder.fill(9, (float)mcMom.y());
-	    hdmcClusFolder.fill(10, (float)mcMom.z());
-	    hdmcClusFolder.fill(11,(float) vert.x());
-	    hdmcClusFolder.fill(12,(float) vert.y());
-	    hdmcClusFolder.fill(13,(float) vert.z());
-	    hdmcClusFolder.addRow();
+	  MCParticle mcpart = getMCParticleInCluster(clust);
+	  if(mcpart!=null) {
+	    float charge = (float)mcpart.getCharge();
+	    int particleID = mcpart.getType().getPDGID();
+	    float mcE = (float)mcpart.getEnergy();
+	    Hep3Vector mcMom = mcpart.getMomentum();
+	    double[] pos = clust.getPosition();
+	    Hep3Vector vert = mcpart.getOrigin();
+
+	    ITuple mcClusFolder = null;
+	    if(_calType=="EM") mcClusFolder = emmcClusFolder;
+	    if(_calType=="HD") mcClusFolder = hdmcClusFolder;
+	    if(mcClusFolder!=null) {
+	      mcClusFolder.fill(0,(float) clust.getRawEnergy());
+	      mcClusFolder.fill(1, (float)pos[0] );
+	      mcClusFolder.fill(2, (float)pos[1] );
+	      mcClusFolder.fill(3, (float)pos[2] );
+	      mcClusFolder.fill(4, clust.getSize() );
+	      mcClusFolder.fill(5, particleID);
+	      mcClusFolder.fill(6, charge);
+	      mcClusFolder.fill(7, mcE);
+	      mcClusFolder.fill(8, (float)mcMom.x());
+	      mcClusFolder.fill(9,(float)mcMom.y());
+	      mcClusFolder.fill(10,(float)mcMom.z());
+	      mcClusFolder.fill(11,(float) vert.x());
+	      mcClusFolder.fill(12,(float) vert.y());
+	      mcClusFolder.fill(13,(float) vert.z());
+	      mcClusFolder.addRow();
+	    }
+	    tag++;
 	  }
-	  tag++;
-        }
+	}
 
         if(_calType=="EM"){
-          _tupleEM.fill(0,recon.size());
+	  int nreco = recon.size();
+          _tupleEM.fill(0,nreco);
 	  _tupleEM.fill(3,_evtnum);
 	  _tupleEM.addRow();
 	  _tupleMCem.fill(0,mctruth.size());
 	  _tupleMCem.fill(3,_evtnum);
 	  _tupleMCem.addRow();
+	  System.out.println("Adding nclusEM="+recon.size());
 	  _tupleEvt.fill(0,_evtnum);
 	  _tupleEvt.addRow();
         }
@@ -316,15 +304,26 @@
 	  _tupleMChd.fill(0,mctruth.size());
 	  _tupleMChd.fill(3,_evtnum);
 	  _tupleMChd.addRow();
+	  System.out.println("Adding nclusHD="+recon.size());
         }
 
 	recTagMap.clear();
     }
 
-    private void processCollection(EventHeader event, String colName) {
-
-      Map<Long,CalorimeterHit> hitmap = _loader.getCollHitMap(colName);
+    private void processCollection(EventHeader event, String colName,
+				   Map<Long,CalorimeterHit> hitmap)
+    {
       if(hitmap.size()==0) return;
+      IDDecoder decoder = _expert.getIDDecoder(colName);
+
+      // get NearestNeighbor clusters
+      List<BasicCluster> nnClusColl = new ArrayList<BasicCluster>();
+      try {
+	nnClusColl = event.get(BasicCluster.class, colName+"NNClusters");
+      }
+      catch(IllegalArgumentException x) {
+	  // ignore
+      }
 
       // get reco clusters
       List<BasicCluster> recoClusColl = new ArrayList<BasicCluster>();
@@ -340,8 +339,11 @@
 	  = event.get(BasicCluster.class, colName+"Clusters");
 
       Collections.sort( mcClusColl, new SortClustersBySize() );
-      System.out.println("*** DTreeAnalysis: colName="+colName
-			 +", # of MCclusters: "+mcClusColl.size());
+      if(_debug>0) {
+	System.out.println("*** DTreeAnalysis: colName="+colName
+			   +", # of MCclusters: "+mcClusColl.size()
+			   +", # of NNclusters: "+nnClusColl.size());
+      }
 
       // Analyze MC clusters
       int mchits=0;
@@ -352,11 +354,25 @@
       }
       int nclusCheat = mcClusColl==null ? 0 : mcClusColl.size();
 
-      System.out.println(" Comparing cheat clusters: "
-			 +", cheater="+nclusCheat);
-      System.out.println(" Comparing total #hits: "
-			 +" in event: "+hitmap.size()
-			 +", in MCClusters="+mchits );
+      // Analyze NN clusters
+      int nnhits=0;
+      for( BasicCluster iclus : nnClusColl ) {
+	double pos[] = iclus.getPosition();
+	List<CalorimeterHit> iHits = iclus.getCalorimeterHits();
+	nnhits += iHits.size(); // includes hits below threshold
+      }
+      int nNNclus = nnClusColl==null ? 0 : nnClusColl.size();
+
+      if(_debug>0) {
+	System.out.println(" Comparing cheat clusters: "
+			   +" cheater="+nclusCheat
+			   +" NN="+nNNclus
+			   );
+	System.out.println(" Comparing total #hits: "
+			   +" in event: "+hitmap.size()
+			   +", in MCClusters="+mchits
+			   +", in NNClusters="+nnhits );
+      }
 
       int nhits = hitmap.size();
       AIDA aida = AIDA.defaultInstance();
@@ -365,7 +381,7 @@
       aida.cloud1D(colName+"-diffHitsColl-mc").fill(nhits-mchits);
       aida.cloud1D(colName+"-numCheatClusters").fill(nclusCheat);
 
-      fillDirectedTreeNtuple( recoClusColl, mcClusColl );
+      fillDirectedTreeNtuple( recoClusColl, mcClusColl, decoder );
     }
 
     private MCParticle getMCParticleInCluster(BasicCluster clust) {
@@ -381,6 +397,7 @@
     private static Map<String,Readout> _roMap;
     private static RunControlParameters _runPar;
     private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+    private CalHitMapMgr _expert = CalHitMapMgr.getInstance();
     private Map<Long,Long> _ParentMap;
     private Map<Long,CalorimeterHit> embhitmap,emehitmap,hdbhitmap,hdehitmap;
     private String _distType,_calType;
@@ -400,4 +417,6 @@
     private ITuple emmcClusFolder,emmcCellFolder;
     private ITuple hdmcClusFolder,hdmcCellFolder;
     private ClusterBuilder clusBuilder = new ClusterBuilder();
+    private int layerIndex, uIndex, vIndex;
+    private double _emEcut, _hdEcut;
 }

lcsim/src/org/lcsim/recon/cluster/density
DirectedTreeClusterer.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- DirectedTreeClusterer.java	12 Dec 2005 05:24:46 -0000	1.2
+++ DirectedTreeClusterer.java	27 Dec 2005 16:33:50 -0000	1.3
@@ -14,20 +14,26 @@
 import org.lcsim.geometry.compact.Readout;
 import org.lcsim.geometry.segmentation.SegmentationBase;
 import org.lcsim.digisim.CellSelector;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
 import org.lcsim.recon.cluster.util.BasicCluster;
 import org.lcsim.recon.cluster.util.SortClustersBySize;
 
 public class DirectedTreeClusterer extends Driver {
 
-    public DirectedTreeClusterer() throws IOException {
+    private int _debug = 0;
+
+    public DirectedTreeClusterer() {
         _runPar = RunControlParameters.getInstance();
+	_emEcut = _runPar.getEMmip() * _runPar.getEMthresh() / _runPar.getEMweight();
+	_hdEcut = _runPar.getHDmip() * _runPar.getHDthresh() / _runPar.getHDweight();
     }
 
     public void process(EventHeader event) {
 
-	System.out.println("*************** New Event **********");
         _evtnum = event.getEventNumber();
-	System.out.println("Starting DirectedTreeClusterer for event no. "+_evtnum);
+	if(_debug>0) {
+	  System.out.println("Start DirectedTreeClusterer, event #"+_evtnum);
+	}
 
         _roMap = event.getDetector().getReadouts();
 
@@ -35,24 +41,30 @@
 	String emeName = "EcalEndcapHits";
 	String hdbName = "HcalBarrHits";
 	String hdeName = "HcalEndcapHits";
-	embhitmap = _loader.getCollHitMap(embName);
-	emehitmap = _loader.getCollHitMap(emeName);
-	hdbhitmap = _loader.getCollHitMap(hdbName);
-	hdehitmap = _loader.getCollHitMap(hdeName);
-
-	System.out.println("DTree: #hits: EMB="+embhitmap.size()
-			   +", EMEC="+emehitmap.size()
-			   +", HB="+hdbhitmap.size()
-			   +", HEC="+hdehitmap.size());
+	embhitmap = _expert.getCollHitMap(embName, _emEcut);
+	emehitmap = _expert.getCollHitMap(emeName, _emEcut);
+	hdbhitmap = _expert.getCollHitMap(hdbName, _hdEcut);
+	hdehitmap = _expert.getCollHitMap(hdeName, _hdEcut);
+	_loader.setDensities(embName, embhitmap);
+	_loader.setDensities(emeName, emehitmap);
+	_loader.setDensities(hdbName, hdbhitmap);
+	_loader.setDensities(hdeName, hdehitmap);
+
+	if(_debug>0) {
+	    System.out.println("DTree: #hits: EMB="+embhitmap.size()
+			       +", EMEC="+emehitmap.size()
+			       +", HB="+hdbhitmap.size()
+			       +", HEC="+hdehitmap.size());
+	}
 
-        if(_runPar.ClusterSeparately()){
+        if(_runPar.ClusterSeparately()) {
             _calType = "EM";
-	    processCollection(event,embName);
-// 	    processCollection(event,emeName);
+	    processCollection(event, embName, embhitmap);
+//  	    processCollection(event, emeName, emehitmap);
 
             _calType = "HD";
-	    processCollection(event,hdbName);
-// 	    processCollection(event,hdeName);
+	    processCollection(event, hdbName, hdbhitmap);
+//  	    processCollection(event, hdeName, hdehitmap);
         }
         else {
 	  assert false : "Sorry, single-pass clustering unavailable for now.";
@@ -79,7 +91,9 @@
     public Vector<Vector<CalorimeterHit>>
 	makeTree(Map<Long,CalorimeterHit> cahitmap, SegmentationBase segm) {
 
-	System.out.println("makeTree called: cahitmap size="+cahitmap.size());
+	if(_debug>0) {
+	  System.out.println("makeTree: cahitmap size="+cahitmap.size());
+	}
         _ParentMap = new HashMap<Long,Long>();
         Vector<CalorimeterHit> RootVec = new Vector<CalorimeterHit>();
 
@@ -161,8 +175,6 @@
 //  		System.out.println("Processing ihit: cellID="+MyTools.printID(cellid));
 		long neigh[] = segm.getNeighbourIDs(_nLyr, _nZ, _nPhi);
 
-		//		if(_evtnum==81 && cellid== 1672347904) System.out.println("neighborhood="+" "+_nLyr+" "+_nZ+" "+_nPhi);
-
                 for(int j=0; j<neigh.length; j++){
                     long jid = neigh[j];
                     CalorimeterHit jhit = cahitmap.get(jid);
@@ -229,13 +241,10 @@
                         double dmin = 9999.;
                         long parentKey = 0;
                         for(Long jkey : nVec) {
-//			    if(_evtnum==81 && cellid== 1672347904) System.out.println(jkey);
                             CalorimeterHit jhit = cahitmap.get(jkey);
                             double[] jpos = jhit.getPosition();
                             CalculateDistance calcD = new CalculateDistance();
                             double distance = calcD.CalculateDistance(_distType,ipos,jpos);
-//			    if(_evtnum==81 && cellid== 1672347904) System.out.println(jkey+" "+myhit2+" "+distance+" "+
-//			       ipos[0]+" "+ipos[1]+" "+ipos[2]+" "+jpos[0]+" "+jpos[1]+" "+jpos[2]);
                             if(distance<dmin){
                                 dmin = distance;
                                 parentKey = jkey;
@@ -261,7 +270,9 @@
                 startingPoints.add(ihit);
             }
         }
-	System.out.println("# starting points = "+startingPoints.size());
+	if(_debug>0) {
+	    System.out.println("# starting points = "+startingPoints.size());
+	}
 
         Vector<Vector<CalorimeterHit>> branches = new Vector<Vector<CalorimeterHit>>();
         for(int i=0;i<startingPoints.size();i++) {
@@ -282,7 +293,7 @@
 
 	int nTrivial = 0;
         int vsiz = RootVec.size();
-	System.out.println("no. of roots = "+vsiz);
+	if(_debug>0) System.out.println("no. of roots = "+vsiz);
         Vector<Vector<CalorimeterHit>> caTrees = new Vector<Vector<CalorimeterHit>>();
         for(int i=0;i<vsiz;i++){
             caTrees.add( new Vector<CalorimeterHit>() );
@@ -305,13 +316,17 @@
 	      }
 	    }
 	  }
-	  if(caTrees.get(i).size()>1) {
-	      System.out.println("ROOT id="+MyTools.printID(tmpID)
-				 +", #hits="+caTrees.get(i).size());
+	  if(_debug>0) {
+	      if(caTrees.get(i).size()>1) {
+		  System.out.println("ROOT id="+MyTools.printID(tmpID)
+				     +", #hits="+caTrees.get(i).size());
+	      }
+	      else 	      ++nTrivial;  // for single-hit clusters
 	  }
-	  else 	      ++nTrivial;  // for single-hit clusters
         }
-	System.out.println("*** plus "+nTrivial+" single-hit clusters.");
+	if(_debug>0) {
+	  System.out.println("*** plus "+nTrivial+" single-hit clusters.");
+	}
 
 // 	for( Vector<MyCalorimeterHit> rootHits : caTrees ) {
 // 	    System.out.println("ROOT: #hits="+rootHits.size());
@@ -323,23 +338,28 @@
     public int emDensityRegion(double dens,int nl,int nz,int np){
 	int region = 0;
 	//	if(nl==6 && nz==3 && np==3){
-	if(dens<5.) region=1;
-	if(dens>=5. && dens<25.) region=2;
-	if(dens>=25.) region=3;
+	if(dens<6.) region=1;
+	if(dens>=6. && dens<26.) region=2;
+	if(dens>=26.) region=3;
 	//	}
 	return region;
     }
 
     public int hdDensityRegion(double dens,int nl,int nz,int np){
 	int region = 0;
+	if(dens<=6.) region=1;
+	if(dens>=6. && dens<26.) region=2;
+	if(dens>=26.) region=3;
 	return region;
     }
 
-    private void processCollection(EventHeader event, String colName) {
+    private void processCollection(EventHeader event, String colName,
+				   Map<Long,CalorimeterHit> hitmap) {
 
-      Map<Long,CalorimeterHit> hitmap = _loader.getCollHitMap(colName);
-      System.out.println("*** procColl: colName="+colName
-			 +", # hits="+hitmap.size());
+      if(_debug>0) {
+	  System.out.println("*** procColl: colName="+colName
+			     +", # hits="+hitmap.size());
+      }
 
       SegmentationBase segm
 	  = (SegmentationBase)_roMap.get(colName).getSegmentation();
@@ -350,7 +370,6 @@
 	  recoClusColl = clusBuilder.makeClusters( trees );
       if (recoClusColl.size() > 0) {
 	  String newName = new String(colName+"DTreeClusters");
-	  System.out.println("Appending to event: <"+newName+">");
 	  event.put( newName, recoClusColl );
       }
 
@@ -365,6 +384,7 @@
     private static Map<String,Readout> _roMap;
     private static RunControlParameters _runPar;
     private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+    private CalHitMapMgr _expert = CalHitMapMgr.getInstance();
     private Map<Long,Long> _ParentMap;
     private Map<Long,CalorimeterHit> embhitmap,emehitmap,hdbhitmap,hdehitmap;
     private String _distType,_calType;
@@ -372,6 +392,7 @@
     private int _nZ;
     private int _nPhi;
     private int _evtnum;
+    private double _emEcut, _hdEcut;
 
     private ClusterBuilder clusBuilder = new ClusterBuilder();
 }
CVSspam 0.2.8