2 added + 1 removed, total 3 files
lcsim/src/org/lcsim/recon/cluster/density
diff -N DTreeAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DTreeAnalysis.java 7 Dec 2005 18:57:14 -0000 1.1
@@ -0,0 +1,364 @@
+package org.lcsim.recon.cluster.density;
+
+import java.io.IOException;
+import java.util.*;
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimCalorimeterHit;
+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.BasicCluster;
+import org.lcsim.recon.cluster.util.SortClustersBySize;
+
+public class DTreeAnalysis extends Driver {
+
+ public DTreeAnalysis() throws IOException {
+ _runPar = RunControlParameters.getInstance();
+ AIDA aida = AIDA.defaultInstance();
+ _tree = aida.tree();
+ _tf = aida.analysisFactory().createTupleFactory(_tree);
+ _tree.mkdir("DirectedTree");
+ _tree.cd("DirectedTree");
+ bookDirectedTreeNtuple();
+ _tree.cd("..");
+ }
+
+ public void process(EventHeader event) {
+
+ System.out.println("*************** New Event **********");
+ _evtnum = event.getEventNumber();
+ System.out.println("Starting DTreeAnalysis for event no. "+_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());
+
+ if(_runPar.ClusterSeparately()){
+ _calType = "EM";
+ processCollection(event,embName);
+// processCollection(event,emeName);
+
+ _calType = "HD";
+ processCollection(event,hdbName);
+// processCollection(event,hdeName);
+ }
+ else {
+ assert false : "Sorry, single-pass clustering unavailable for now.";
+ }
+
+ _tree.cd("..");
+ }
+
+ public void bookDirectedTreeNtuple(){
+
+ String[] columnNamesEvent = {"evtno"};
+ Class[] columnClassesEvent = {Integer.TYPE};
+ _tupleEvt = _tf.create("EVT","Event Parameters",columnNamesEvent,columnClassesEvent);
+
+ if(_runPar.ClusterSeparately()){
+ String[] columnNamesEM = {"nclusEM","emClusFolder={float ene,float x,float y,float z,int size}",
+ "emCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTag}","evtno"};
+ Class[] columnClassesEM = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
+ _tupleEM = _tf.create("EM","Recon Clusters",columnNamesEM,columnClassesEM);
+
+ String[] columnNamesHD = {"nclusHD","hdClusFolder={float ene,float x,float y,float z,int size}",
+ "hdCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTag}","evtno"};
+ Class[] columnClassesHD = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
+ _tupleHD = _tf.create("HD","Recon Clusters",columnNamesHD,columnClassesHD);
+
+ String[] columnNamesMCem = {"nclusMCem","emmcClusFolder={float ene,float x,float y,float z,int size,int pid,float chrg,float e,float px,float py,float pz}","emmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi}","evtno"};
+ 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}","hdmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi}","evtno"};
+ Class[] columnClassesMChd = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
+ _tupleMChd = _tf.create("MCHD","Gen Clusters in HD",columnNamesMChd,columnClassesMChd);
+ }
+ else{
+ }
+ }
+
+ public void fillDirectedTreeNtuple(List<BasicCluster> recon,
+ List<BasicCluster> mctruth)
+ {
+ if(_calType=="EM"){
+ emClusFolder = _tupleEM.getTuple(1);
+ emCellFolder = _tupleEM.getTuple(2);
+ emmcClusFolder = _tupleMCem.getTuple(1);
+ emmcCellFolder = _tupleMCem.getTuple(2);
+ }
+
+ if(_calType=="HD"){
+ hdClusFolder = _tupleHD.getTuple(1);
+ hdCellFolder = _tupleHD.getTuple(2);
+ hdmcClusFolder = _tupleMChd.getTuple(1);
+ hdmcCellFolder = _tupleMChd.getTuple(2);
+ }
+
+ 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("Damn 3");
+// for(int i=0;i<cellVec.size();i++){
+// MyCalorimeterHit myhit = (MyCalorimeterHit) cellVec.get(i);
+ for(CalorimeterHit ihit : cellVec) {
+ Integer tagObject = new Integer(tag);
+// System.out.println("putting <"+myhit+"->"+tagObject+">");
+ recTagMap.put(ihit,tagObject);
+ double[] pos = ihit.getPosition();
+ double cdens = _loader.getDensity(ihit);
+ 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.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.addRow();
+ }
+ }
+ ITuple clusFolder = null;
+ if(_calType=="EM") clusFolder = emClusFolder;
+ if(_calType=="HD") clusFolder = hdClusFolder;
+ if(clusFolder!=null) {
+ Hep3Vector pos = new BasicHep3Vector( clust.getPosition() );
+ clusFolder.fill(0,(float) clust.getEnergy());
+ 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++;
+ }
+
+ tag = 0;
+ for( BasicCluster clust : mctruth ) {
+ for( CalorimeterHit ihit : clust.getCalorimeterHits() ) {
+ Integer tagObject = recTagMap.get(ihit);
+// System.out.println("getting: "+ihit+" -> "+tagObject);
+
+ int rtag = -1;
+ if(tagObject!=null) rtag = tagObject.intValue();
+ double[] pos = ihit.getPosition();
+ double cdens = _loader.getDensity(ihit);
+ if(cdens==0) continue; // skip hits below energy cut
+ long jid = ihit.getCellID();
+ int ily = MyTools.getLayer(jid);
+ int iz = MyTools.getThetaBin(jid);
+ int iphi = MyTools.getPhiBin(jid);
+ 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.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.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();
+ 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.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.addRow();
+ }
+ tag++;
+ }
+
+ if(_calType=="EM"){
+ _tupleEM.fill(0,recon.size());
+ _tupleEM.fill(3,_evtnum);
+ _tupleEM.addRow();
+ _tupleMCem.fill(0,mctruth.size());
+ _tupleMCem.fill(3,_evtnum);
+ _tupleMCem.addRow();
+ _tupleEvt.fill(0,_evtnum);
+ _tupleEvt.addRow();
+ }
+ if(_calType=="HD"){
+ _tupleHD.fill(0,recon.size());
+ _tupleHD.fill(3,_evtnum);
+ _tupleHD.addRow();
+ _tupleMChd.fill(0,mctruth.size());
+ _tupleMChd.fill(3,_evtnum);
+ _tupleMChd.addRow();
+ }
+
+ recTagMap.clear();
+ }
+
+ private void processCollection(EventHeader event, String colName) {
+
+ Map<Long,CalorimeterHit> hitmap = _loader.getCollHitMap(colName);
+ if(hitmap.size()==0) return;
+
+ // get reco clusters
+ List<BasicCluster> recoClusColl = new ArrayList<BasicCluster>();
+ try {
+ recoClusColl = event.get(BasicCluster.class, colName+"DTreeClusters");
+ }
+ catch(IllegalArgumentException x) {
+ // ignore
+ }
+
+ // get MC clusters
+ List<BasicCluster> mcClusColl
+ = event.get(BasicCluster.class, colName+"Clusters");
+
+ Collections.sort( mcClusColl, new SortClustersBySize() );
+ System.out.println("*** DTreeAnalysis: colName="+colName
+ +", # of MCclusters: "+mcClusColl.size());
+
+ // Analyze MC clusters
+ int mchits=0;
+ for( BasicCluster iclus : mcClusColl ) {
+ double pos[] = iclus.getPosition();
+ List<CalorimeterHit> iHits = iclus.getCalorimeterHits();
+ mchits += iHits.size(); // includes hits below threshold
+ }
+ 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 );
+
+ int nhits = hitmap.size();
+ AIDA aida = AIDA.defaultInstance();
+ aida.cloud1D(colName+"-Nhits").fill( nhits );
+ aida.cloud1D(colName+"-NmcHits").fill(mchits);
+ aida.cloud1D(colName+"-diffHitsColl-mc").fill(nhits-mchits);
+ aida.cloud1D(colName+"-numCheatClusters").fill(nclusCheat);
+
+ fillDirectedTreeNtuple( recoClusColl, mcClusColl );
+ }
+
+ private MCParticle getMCParticleInCluster(BasicCluster clust) {
+ for( CalorimeterHit hit : clust.getCalorimeterHits() ) {
+ SimCalorimeterHit simhit = (SimCalorimeterHit)hit;
+ if( simhit!=null && simhit.getMCParticleCount()==1 ) {
+ return simhit.getMCParticle(0);
+ }
+ }
+ return null;
+ }
+
+ private static Map<String,Readout> _roMap;
+ private static RunControlParameters _runPar;
+ private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+ private Map<Long,Long> _ParentMap;
+ private Map<Long,CalorimeterHit> embhitmap,emehitmap,hdbhitmap,hdehitmap;
+ private String _distType,_calType;
+ private int _nLyr;
+ private int _nZ;
+ private int _nPhi;
+ private int _evtnum;
+
+ private ITree _tree;
+ private ITupleFactory _tf;
+ private ITuple _tupleEM;
+ private ITuple _tupleHD;
+ private ITuple _tupleMCem,_tupleMChd;
+ private ITuple _tupleEvt;
+ private ITuple emClusFolder,emCellFolder;
+ private ITuple hdClusFolder,hdCellFolder;
+ private ITuple emmcClusFolder,emmcCellFolder;
+ private ITuple hdmcClusFolder,hdmcCellFolder;
+ private ClusterBuilder clusBuilder = new ClusterBuilder();
+}
lcsim/src/org/lcsim/recon/cluster/density
diff -N DirectedTreeClusterer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DirectedTreeClusterer.java 7 Dec 2005 18:57:14 -0000 1.1
@@ -0,0 +1,378 @@
+package org.lcsim.recon.cluster.density;
+
+import java.io.IOException;
+import java.util.*;
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+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.BasicCluster;
+import org.lcsim.recon.cluster.util.SortClustersBySize;
+
+public class DirectedTreeClusterer extends Driver {
+
+ public DirectedTreeClusterer() throws IOException {
+ _runPar = RunControlParameters.getInstance();
+ }
+
+ public void process(EventHeader event) {
+
+ System.out.println("*************** New Event **********");
+ _evtnum = event.getEventNumber();
+ System.out.println("Starting DirectedTreeClusterer for event no. "+_evtnum);
+
+ _roMap = event.getDetector().getReadouts();
+
+ 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());
+
+ if(_runPar.ClusterSeparately()){
+ _calType = "EM";
+ processCollection(event,embName);
+// processCollection(event,emeName);
+
+ _calType = "HD";
+ processCollection(event,hdbName);
+// processCollection(event,hdeName);
+ }
+ else {
+ assert false : "Sorry, single-pass clustering unavailable for now.";
+// Vector calhit = new Vector();
+// calhit.addAll(embhit.values());
+// calhit.addAll(emehit.values());
+// calhit.addAll(hdbhit.values());
+// calhit.addAll(hdehit.values());
+
+// cahitmap = new HashMap();
+// cahitmap.putAll(embhitmap);
+// cahitmap.putAll(emehitmap);
+// cahitmap.putAll(hdbhitmap);
+// cahitmap.putAll(hdehitmap);
+
+// Vector[] calTrees = makeTree(calhit,cahitmap);
+
+// ClusterBuilder clus = new ClusterBuilder();
+// List caClusters = clus.makeClusters(calTrees,_runPar);
+ }
+
+ }
+
+ public Vector<Vector<CalorimeterHit>>
+ makeTree(Map<Long,CalorimeterHit> cahitmap, SegmentationBase segm) {
+
+ System.out.println("makeTree called: cahitmap size="+cahitmap.size());
+ _ParentMap = new HashMap<Long,Long>();
+ Vector<CalorimeterHit> RootVec = new Vector<CalorimeterHit>();
+
+ for( CalorimeterHit ihit : cahitmap.values() ) {
+ long cellid = ihit.getCellID();
+ double idens = _loader.getDensity( cellid );
+// if(_evtnum==81 && cellid== 1672347904) System.out.println("density="+idens);
+
+ segm.setID(cellid);
+// System.out.println("makeTree ihit: cellID="+MyTools.printID(cellid)
+// +", layer "+segm.getLayer()
+// +", dens="+idens);
+
+ double maxdensDiff = -999999.;
+ long maxdensID = -999999;
+ Vector<Long> nVec = new Vector<Long>();
+ if(idens==0.0) {
+ RootVec.add(ihit);
+// System.out.println("Zero density root");
+ }
+ else{
+ double[] ipos = ihit.getPosition();
+ if(_calType=="EM"){
+ int nLyrOrig = _runPar.getLyrNeighEM();
+ int nZOrig = _runPar.getZNeighEM();
+ int nPhiOrig = _runPar.getPhiNeighEM();
+ int dRegion = emDensityRegion(idens,nLyrOrig,nZOrig,nPhiOrig);
+ int[] lyrCon = _runPar.getLyrContracEM();
+ int[] zCon = _runPar.getZContracEM();
+ int[] phiCon = _runPar.getPhiContracEM();
+// System.out.println("Original Window"+" "+nLyrOrig+" "+nZOrig+" "+nPhiOrig);
+ if(dRegion==1){
+ _nLyr = nLyrOrig - lyrCon[0];
+ _nZ = nZOrig - zCon[0];
+ _nPhi = nPhiOrig - phiCon[0];
+ }
+ if(dRegion==2){
+ _nLyr = nLyrOrig - lyrCon[1];
+ _nZ = nZOrig - zCon[1];
+ _nPhi = nPhiOrig - phiCon[1];
+ }
+ if(dRegion==3){
+ _nLyr = nLyrOrig - lyrCon[2];
+ _nZ = nZOrig - zCon[2];
+ _nPhi = nPhiOrig - phiCon[2];
+ }
+ }
+ if(_calType=="HD"){
+ int nLyrOrig = _runPar.getLyrNeighHD();
+ int nZOrig = _runPar.getZNeighHD();
+ int nPhiOrig = _runPar.getPhiNeighHD();
+ int dRegion = hdDensityRegion(idens,nLyrOrig,nZOrig,nPhiOrig);
+ int[] lyrCon = _runPar.getLyrContracHD();
+ int[] zCon = _runPar.getZContracHD();
+ int[] phiCon = _runPar.getPhiContracHD();
+// System.out.println("Original Window"+" "+nLyrOrig+" "+nZOrig+" "+nPhiOrig);
+ if(dRegion==0){
+ _nLyr = nLyrOrig;
+ _nZ = nZOrig;
+ _nPhi = nPhiOrig;
+ }
+ if(dRegion==1){
+ _nLyr = nLyrOrig - lyrCon[0];
+ _nZ = nZOrig - zCon[0];
+ _nPhi = nPhiOrig - phiCon[0];
+ }
+ if(dRegion==2){
+ _nLyr = nLyrOrig - lyrCon[1];
+ _nZ = nZOrig - zCon[1];
+ _nPhi = nPhiOrig - phiCon[1];
+ }
+ if(dRegion==3){
+ _nLyr = nLyrOrig - lyrCon[2];
+ _nZ = nZOrig - zCon[2];
+ _nPhi = nPhiOrig - phiCon[2];
+ }
+ }
+
+// int neigh[] = segm.getNeighbors(cellid,_nLyr,_nZ,_nPhi);
+// if(_calType=="HD")
+// System.out.println("window in directed tree "+_nLyr+" "+_nZ+" "+_nPhi);
+ segm.setID(cellid);
+// 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);
+ if(jhit==null) continue;
+ double jdens = _loader.getDensity(jhit);
+ double[] jpos = jhit.getPosition();
+
+ _distType = _runPar.getDistanceType();
+ CalculateDistance calcD = new CalculateDistance();
+ double distance = calcD.CalculateDistance(_distType,ipos,jpos);
+ double densDiff = (jdens-idens)/distance;
+//??? if(distance==0) densDiff=0;
+
+ if(densDiff==0.0) nVec.add(jid);
+// System.out.println(" jhit: layer="+(jid&0x7f)
+// +", cellID="+MyTools.printID(jid)
+// +", idens="+idens
+// +", jdens="+jdens
+// +", dist="+distance
+// +", densDiff="+densDiff);
+
+ if(densDiff>maxdensDiff){
+ maxdensDiff = densDiff;
+ maxdensID = jid;
+ }
+ calcD = null;
+ }
+
+// if(maxdensDiff<0.0){
+// RootVec.add(myhit1);
+// System.out.println("Bonafide root");
+// }
+
+ if(maxdensDiff>0.0){
+ Long key1 = new Long(cellid);
+ Long key2 = new Long(maxdensID);
+ _ParentMap.put(key1,key2);
+// System.out.println(" maxdensDiff>0: "
+// +" id="+MyTools.printID(maxdensID)
+// +", parent of "+MyTools.printID(cellid));
+ }
+
+ if(maxdensDiff==0.0){
+ Vector<Long> removeItems = new Vector<Long>();
+// Long key = new Long(cellid);
+ for(Long jkey : nVec) {
+ Vector<Long> temporary = new Vector<Long>();
+ Long parent = _ParentMap.get(jkey);
+
+ while(parent!=null){
+ temporary.add(parent);
+ Long lkey = parent;
+ parent = _ParentMap.get(lkey);
+ }
+ if(temporary.contains(cellid)){
+ removeItems.add(jkey);
+ }
+ }
+ nVec.removeAll(removeItems);
+
+ if(nVec.size()==0){
+ RootVec.add(ihit);
+ }
+ else{
+ 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;
+ }
+ }
+ _ParentMap.put(cellid,parentKey);
+// System.out.println("Cell is matched to identical density neighbor");
+ }
+ }
+ }
+ }
+
+ Vector<CalorimeterHit> startingPoints = new Vector<CalorimeterHit>();
+ for( CalorimeterHit ihit : cahitmap.values() ) {
+ long cellid = ihit.getCellID();
+ boolean isParent = _ParentMap.containsValue(cellid);
+ boolean isRoot = RootVec.contains(ihit);
+ if(!isParent && !isRoot){
+ startingPoints.add(ihit);
+ }
+ }
+ System.out.println("# starting points = "+startingPoints.size());
+
+ Vector<Vector<CalorimeterHit>> branches = new Vector<Vector<CalorimeterHit>>();
+ for(int i=0;i<startingPoints.size();i++) {
+ branches.add( new Vector<CalorimeterHit>() );
+ }
+ for(int i=0;i<startingPoints.size();i++){
+ CalorimeterHit ihit = startingPoints.get(i);
+ branches.get(i).add(ihit);
+ long cellid = ihit.getCellID();
+ Long key = new Long(cellid);
+ while(_ParentMap.containsKey(cellid)){
+ Long jkey = _ParentMap.get(cellid);
+ CalorimeterHit jhit = cahitmap.get(jkey);
+ branches.get(i).add(jhit);
+ cellid = jkey.longValue();
+ }
+ }
+
+ int nTrivial = 0;
+ int vsiz = RootVec.size();
+ 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>() );
+ }
+
+ for(int i=0;i<vsiz;i++){
+ CalorimeterHit ihit = RootVec.get(i);
+ long tmpID = ihit.getCellID();
+ caTrees.get(i).add(ihit);
+ for(int j=0;j<startingPoints.size();j++) {
+ CalorimeterHit parentHit = branches.get(j).lastElement();
+// if(_evtnum==81) System.out.println(myhit1+" "+parentHit+" "+startingPoints.size()+" "+j);
+ if(parentHit.equals(ihit)) {
+ for( CalorimeterHit jhit : branches.get(j) ) {
+ if(!caTrees.get(i).contains(jhit)){
+ caTrees.get(i).add(jhit);
+ long tmppID = jhit.getCellID();
+// System.out.println("Branches id="+MyTools.printID(tmppID));
+ }
+ }
+ }
+ }
+ 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
+ }
+ System.out.println("*** plus "+nTrivial+" single-hit clusters.");
+
+// for( Vector<MyCalorimeterHit> rootHits : caTrees ) {
+// System.out.println("ROOT: #hits="+rootHits.size());
+// }
+
+ return caTrees;
+ }
+
+ 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;
+ // }
+ return region;
+ }
+
+ public int hdDensityRegion(double dens,int nl,int nz,int np){
+ int region = 0;
+ return region;
+ }
+
+ private void processCollection(EventHeader event, String colName) {
+
+ Map<Long,CalorimeterHit> hitmap = _loader.getCollHitMap(colName);
+ System.out.println("*** procColl: colName="+colName
+ +", # hits="+hitmap.size());
+
+ SegmentationBase segm
+ = (SegmentationBase)_roMap.get(colName).getSegmentation();
+ Vector<Vector<CalorimeterHit>> trees = makeTree(hitmap, segm);
+
+ List<BasicCluster> recoClusColl = new ArrayList<BasicCluster>();
+ if(trees.size()>0)
+ 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 );
+ }
+
+// int nhits = hitmap.size();
+// AIDA aida = AIDA.defaultInstance();
+// aida.cloud1D(colName+"-Nhits").fill( nhits );
+// aida.cloud1D(colName+"-NmcHits").fill(mchits);
+// aida.cloud1D(colName+"-diffHitsColl-mc").fill(nhits-mchits);
+// aida.cloud1D(colName+"-numCheatClusters").fill(nclusCheat);
+ }
+
+ private static Map<String,Readout> _roMap;
+ private static RunControlParameters _runPar;
+ private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+ private Map<Long,Long> _ParentMap;
+ private Map<Long,CalorimeterHit> embhitmap,emehitmap,hdbhitmap,hdehitmap;
+ private String _distType,_calType;
+ private int _nLyr;
+ private int _nZ;
+ private int _nPhi;
+ private int _evtnum;
+
+ private ClusterBuilder clusBuilder = new ClusterBuilder();
+}
lcsim/src/org/lcsim/recon/cluster/density
diff -N DirectedTree.java
--- DirectedTree.java 29 Sep 2005 19:08:36 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,632 +0,0 @@
-package org.lcsim.recon.cluster.density;
-
-import java.io.IOException;
-import java.util.*;
-import hep.aida.*;
-import hep.physics.vec.Hep3Vector;
-
-import org.lcsim.util.aida.AIDA;
-import org.lcsim.util.Driver;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.SimCalorimeterHit;
-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.BasicCluster;
-import org.lcsim.recon.cluster.util.SortClustersBySize;
-
-public class DirectedTree extends Driver {
-
- public DirectedTree(ITupleFactory tf,ITree tree) throws IOException{
-
- _runPar = RunControlParameters.getInstance();
-
- _tree = tree;
- _tf = tf;
- _tree.mkdir("DirectedTree");
- _tree.cd("DirectedTree");
- bookDirectedTreeNtuple();
- _tree.cd("..");
- }
-
- public void process(EventHeader event) {
-
- System.out.println("*************** New Event **********");
- AIDA aida = AIDA.defaultInstance();
- _evtnum = event.getEventNumber();
- System.out.println("Starting DirectedTree for event no. "+_evtnum);
- _tree.cd("DirectedTree");
-
- _roMap = event.getDetector().getReadouts();
-
- 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());
-
- if(_runPar.ClusterSeparately()){
- _calType = "EM";
- processCollection(event,embName);
-
- _calType = "HD";
- processCollection(event,hdbName);
- }
- else {
- assert false : "Sorry, single-pass clustering unavailable for now.";
-// Vector calhit = new Vector();
-// calhit.addAll(embhit.values());
-// calhit.addAll(emehit.values());
-// calhit.addAll(hdbhit.values());
-// calhit.addAll(hdehit.values());
-
-// cahitmap = new HashMap();
-// cahitmap.putAll(embhitmap);
-// cahitmap.putAll(emehitmap);
-// cahitmap.putAll(hdbhitmap);
-// cahitmap.putAll(hdehitmap);
-
-// Vector[] calTrees = makeTree(calhit,cahitmap);
-
-// ClusterBuilder clus = new ClusterBuilder();
-// List caClusters = clus.makeClusters(calTrees,_runPar);
- }
-
- _tree.cd("..");
- }
-
- public Vector<Vector<CalorimeterHit>>
- makeTree(Map<Long,CalorimeterHit> cahitmap, SegmentationBase segm) {
-
- System.out.println("makeTree called: cahitmap size="+cahitmap.size());
- _ParentMap = new HashMap<Long,Long>();
- Vector<CalorimeterHit> RootVec = new Vector<CalorimeterHit>();
-
- for( CalorimeterHit ihit : cahitmap.values() ) {
- long cellid = ihit.getCellID();
- double idens = _loader.getDensity( cellid );
-// if(_evtnum==81 && cellid== 1672347904) System.out.println("density="+idens);
-
- segm.setID(cellid);
-// System.out.println("makeTree ihit: cellID="+MyTools.printID(cellid)
-// +", layer "+segm.getLayer()
-// +", dens="+idens);
-
- double maxdensDiff = -999999.;
- long maxdensID = -999999;
- Vector<Long> nVec = new Vector<Long>();
- if(idens==0.0) {
- RootVec.add(ihit);
-// System.out.println("Zero density root");
- }
- else{
- double[] ipos = ihit.getPosition();
- if(_calType=="EM"){
- int nLyrOrig = _runPar.getLyrNeighEM();
- int nZOrig = _runPar.getZNeighEM();
- int nPhiOrig = _runPar.getPhiNeighEM();
- int dRegion = emDensityRegion(idens,nLyrOrig,nZOrig,nPhiOrig);
- int[] lyrCon = _runPar.getLyrContracEM();
- int[] zCon = _runPar.getZContracEM();
- int[] phiCon = _runPar.getPhiContracEM();
-// System.out.println("Original Window"+" "+nLyrOrig+" "+nZOrig+" "+nPhiOrig);
- if(dRegion==1){
- _nLyr = nLyrOrig - lyrCon[0];
- _nZ = nZOrig - zCon[0];
- _nPhi = nPhiOrig - phiCon[0];
- }
- if(dRegion==2){
- _nLyr = nLyrOrig - lyrCon[1];
- _nZ = nZOrig - zCon[1];
- _nPhi = nPhiOrig - phiCon[1];
- }
- if(dRegion==3){
- _nLyr = nLyrOrig - lyrCon[2];
- _nZ = nZOrig - zCon[2];
- _nPhi = nPhiOrig - phiCon[2];
- }
- }
- if(_calType=="HD"){
- int nLyrOrig = _runPar.getLyrNeighHD();
- int nZOrig = _runPar.getZNeighHD();
- int nPhiOrig = _runPar.getPhiNeighHD();
- int dRegion = hdDensityRegion(idens,nLyrOrig,nZOrig,nPhiOrig);
- int[] lyrCon = _runPar.getLyrContracHD();
- int[] zCon = _runPar.getZContracHD();
- int[] phiCon = _runPar.getPhiContracHD();
-// System.out.println("Original Window"+" "+nLyrOrig+" "+nZOrig+" "+nPhiOrig);
- if(dRegion==0){
- _nLyr = nLyrOrig;
- _nZ = nZOrig;
- _nPhi = nPhiOrig;
- }
- if(dRegion==1){
- _nLyr = nLyrOrig - lyrCon[0];
- _nZ = nZOrig - zCon[0];
- _nPhi = nPhiOrig - phiCon[0];
- }
- if(dRegion==2){
- _nLyr = nLyrOrig - lyrCon[1];
- _nZ = nZOrig - zCon[1];
- _nPhi = nPhiOrig - phiCon[1];
- }
- if(dRegion==3){
- _nLyr = nLyrOrig - lyrCon[2];
- _nZ = nZOrig - zCon[2];
- _nPhi = nPhiOrig - phiCon[2];
- }
- }
-
-// int neigh[] = segm.getNeighbors(cellid,_nLyr,_nZ,_nPhi);
-// if(_calType=="HD")
-// System.out.println("window in directed tree "+_nLyr+" "+_nZ+" "+_nPhi);
- segm.setID(cellid);
-// 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);
- if(jhit==null) continue;
- double jdens = _loader.getDensity(jhit);
- double[] jpos = jhit.getPosition();
-
- _distType = _runPar.getDistanceType();
- CalculateDistance calcD = new CalculateDistance();
- double distance = calcD.CalculateDistance(_distType,ipos,jpos);
- double densDiff = (jdens-idens)/distance;
-//??? if(distance==0) densDiff=0;
-
- if(densDiff==0.0) nVec.add(jid);
-// System.out.println(" jhit: layer="+(jid&0x7f)
-// +", cellID="+MyTools.printID(jid)
-// +", idens="+idens
-// +", jdens="+jdens
-// +", dist="+distance
-// +", densDiff="+densDiff);
-
- if(densDiff>maxdensDiff){
- maxdensDiff = densDiff;
- maxdensID = jid;
- }
- calcD = null;
- }
-
-// if(maxdensDiff<0.0){
-// RootVec.add(myhit1);
-// System.out.println("Bonafide root");
-// }
-
- if(maxdensDiff>0.0){
- Long key1 = new Long(cellid);
- Long key2 = new Long(maxdensID);
- _ParentMap.put(key1,key2);
-// System.out.println(" maxdensDiff>0: "
-// +" id="+MyTools.printID(maxdensID)
-// +", parent of "+MyTools.printID(cellid));
- }
-
- if(maxdensDiff==0.0){
- Vector<Long> removeItems = new Vector<Long>();
-// Long key = new Long(cellid);
- for(Long jkey : nVec) {
- Vector<Long> temporary = new Vector<Long>();
- Long parent = _ParentMap.get(jkey);
-
- while(parent!=null){
- temporary.add(parent);
- Long lkey = parent;
- parent = _ParentMap.get(lkey);
- }
- if(temporary.contains(cellid)){
- removeItems.add(jkey);
- }
- }
- nVec.removeAll(removeItems);
-
- if(nVec.size()==0){
- RootVec.add(ihit);
- }
- else{
- double dmin = 9999.;
- long parentKey = 0;
- for(Long jkey : nVec) {
-// if(_evtnum==81 && cellid== 1672347904) System.out.println(jkey);
- CalorimeterHit jhit = cahitmap.get(jkey);
-// SimCalorimeterHit jhit = myhit2.getHit();
- 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;
- }
- }
- _ParentMap.put(cellid,parentKey);
-// System.out.println("Cell is matched to identical density neighbor");
- }
- }
- }
- }
-
- Vector<CalorimeterHit> startingPoints = new Vector<CalorimeterHit>();
- for( CalorimeterHit ihit : cahitmap.values() ) {
- long cellid = ihit.getCellID();
- boolean isParent = _ParentMap.containsValue(cellid);
- boolean isRoot = RootVec.contains(ihit);
- if(!isParent && !isRoot){
- startingPoints.add(ihit);
- }
- }
- System.out.println("# starting points = "+startingPoints.size());
-
- Vector<Vector<CalorimeterHit>> branches = new Vector<Vector<CalorimeterHit>>();
- for(int i=0;i<startingPoints.size();i++) {
- branches.add( new Vector<CalorimeterHit>() );
- }
- for(int i=0;i<startingPoints.size();i++){
- CalorimeterHit ihit = startingPoints.get(i);
- branches.get(i).add(ihit);
- long cellid = ihit.getCellID();
- Long key = new Long(cellid);
- while(_ParentMap.containsKey(cellid)){
- Long jkey = _ParentMap.get(cellid);
- CalorimeterHit jhit = cahitmap.get(jkey);
- branches.get(i).add(jhit);
- cellid = jkey.longValue();
- }
- }
-
- int vsiz = RootVec.size();
- 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>() );
- }
-
- for(int i=0;i<vsiz;i++){
- CalorimeterHit ihit = RootVec.get(i);
- long tmpID = ihit.getCellID();
- caTrees.get(i).add(ihit);
- for(int j=0;j<startingPoints.size();j++) {
- CalorimeterHit parentHit = branches.get(j).lastElement();
-// if(_evtnum==81) System.out.println(myhit1+" "+parentHit+" "+startingPoints.size()+" "+j);
- if(parentHit.equals(ihit)) {
- for( CalorimeterHit jhit : branches.get(j) ) {
- if(!caTrees.get(i).contains(jhit)){
- caTrees.get(i).add(jhit);
- long tmppID = jhit.getCellID();
-// System.out.println("Branches id="+MyTools.printID(tmppID));
- }
- }
- }
- }
- System.out.println("ROOT id="+MyTools.printID(tmpID)
- +", #hits="+caTrees.get(i).size());
- }
-
-// for( Vector<MyCalorimeterHit> rootHits : caTrees ) {
-// System.out.println("ROOT: #hits="+rootHits.size());
-// }
-
- return caTrees;
- }
-
- 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;
- // }
- return region;
- }
-
- public int hdDensityRegion(double dens,int nl,int nz,int np){
-
- int region = 0;
- return region;
- }
-
- public void bookDirectedTreeNtuple(){
-
- String[] columnNamesEvent = {"evtno"};
- Class[] columnClassesEvent = {Integer.TYPE};
- _tupleEvt = _tf.create("EVT","Event Parameters",columnNamesEvent,columnClassesEvent);
-
- if(_runPar.ClusterSeparately()){
- String[] columnNamesEM = {"nclusEM","emClusFolder={float ene,float x,float y,float z,int size}",
- "emCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTag}","evtno"};
- Class[] columnClassesEM = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
- _tupleEM = _tf.create("EM","Recon Clusters",columnNamesEM,columnClassesEM);
-
- String[] columnNamesHD = {"nclusHD","hdClusFolder={float ene,float x,float y,float z,int size}",
- "hdCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTag}","evtno"};
- Class[] columnClassesHD = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
- _tupleHD = _tf.create("HD","Recon Clusters",columnNamesHD,columnClassesHD);
-
- String[] columnNamesMCem = {"nclusMCem","emmcClusFolder={float ene,float x,float y,float z,int size,int pid,float chrg,float e,float px,float py,float pz}","emmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi}","evtno"};
- 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}","hdmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi}","evtno"};
- Class[] columnClassesMChd = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
- _tupleMChd = _tf.create("MCHD","Gen Clusters in HD",columnNamesMChd,columnClassesMChd);
- }
- else{
- }
- }
-
- public void fillDirectedTreeNtuple(List<Cluster> recon,
- List<BasicCluster> mctruth)
- {
- if(_calType=="EM"){
- emClusFolder = _tupleEM.getTuple(1);
- emCellFolder = _tupleEM.getTuple(2);
- emmcClusFolder = _tupleMCem.getTuple(1);
- emmcCellFolder = _tupleMCem.getTuple(2);
- }
-
- if(_calType=="HD"){
- hdClusFolder = _tupleHD.getTuple(1);
- hdCellFolder = _tupleHD.getTuple(2);
- hdmcClusFolder = _tupleMChd.getTuple(1);
- hdmcCellFolder = _tupleMChd.getTuple(2);
- }
-
- int tag = 0;
- Map<CalorimeterHit,Integer> recTagMap
- = new HashMap<CalorimeterHit,Integer>();
- for(Cluster clust : recon ) {
- List<CalorimeterHit> cellVec = clust.getCalorimeterHits();
- // System.out.println(tag+" "+cellVec.size());
- if(cellVec.size()==0.) System.out.println("Damn 3");
-// for(int i=0;i<cellVec.size();i++){
-// MyCalorimeterHit myhit = (MyCalorimeterHit) cellVec.get(i);
- for(CalorimeterHit ihit : cellVec) {
- Integer tagObject = new Integer(tag);
-// System.out.println("putting <"+myhit+"->"+tagObject+">");
- recTagMap.put(ihit,tagObject);
- double[] pos = ihit.getPosition();
- double cdens = _loader.getDensity(ihit);
- 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.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.addRow();
- }
- }
- if(_calType=="EM"){
- emClusFolder.fill(0,(float) clust._ene);
- emClusFolder.fill(1,(float) clust._xpos);
- emClusFolder.fill(2,(float) clust._ypos);
- emClusFolder.fill(3,(float) clust._zpos);
- emClusFolder.fill(4, clust.getSize());
- emClusFolder.addRow();
-// System.out.println("EM Clus "+tag+" "+clust._xpos+" "+clust._ypos+" "+clust._zpos+" "+ clust._size );
- }
- if(_calType=="HD"){
- hdClusFolder.fill(0,(float) clust._ene);
- 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++;
- }
-
- tag = 0;
- for( BasicCluster clust : mctruth ) {
- for( CalorimeterHit ihit : clust.getCalorimeterHits() ) {
- Integer tagObject = recTagMap.get(ihit);
-// System.out.println("getting: "+ihit+" -> "+tagObject);
-
- int rtag = -1;
- if(tagObject!=null) rtag = tagObject.intValue();
- double[] pos = ihit.getPosition();
- double cdens = _loader.getDensity(ihit);
- if(cdens==0) continue; // skip hits below energy cut
- long jid = ihit.getCellID();
- int ily = MyTools.getLayer(jid);
- int iz = MyTools.getThetaBin(jid);
- int iphi = MyTools.getPhiBin(jid);
- 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.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.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();
- 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.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.addRow();
- }
- tag++;
- }
-
- if(_calType=="EM"){
- _tupleEM.fill(0,recon.size());
- _tupleEM.fill(3,_evtnum);
- _tupleEM.addRow();
- _tupleMCem.fill(0,mctruth.size());
- _tupleMCem.fill(3,_evtnum);
- _tupleMCem.addRow();
- _tupleEvt.fill(0,_evtnum);
- _tupleEvt.addRow();
- }
- if(_calType=="HD"){
- _tupleHD.fill(0,recon.size());
- _tupleHD.fill(3,_evtnum);
- _tupleHD.addRow();
- _tupleMChd.fill(0,mctruth.size());
- _tupleMChd.fill(3,_evtnum);
- _tupleMChd.addRow();
- }
-
- recTagMap.clear();
- }
-
- private void processCollection(EventHeader event, String colName) {
-
- System.out.println("*** Entering procColl: colName="+colName);
- Map<Long,CalorimeterHit> hitmap = _loader.getCollHitMap(colName);
-
- SegmentationBase segm
- = (SegmentationBase)_roMap.get(colName).getSegmentation();
- Vector<Vector<CalorimeterHit>> trees = makeTree(hitmap, segm);
-
- List<Cluster> recoClusColl = clusBuilder.makeClusters( trees, _runPar);
-
- int nclusCheat=0;
- List<BasicCluster> mcClusColl = event.get(BasicCluster.class, colName+"Clusters");
- System.out.println("*** procColl: colName="+colName
- +", # of MCclusters: "+mcClusColl.size());
-
- Collections.sort( mcClusColl, new SortClustersBySize() );
-
- // Analyze MC clusters
- int mchits=0;
- for( BasicCluster iclus : mcClusColl ) {
- double pos[] = iclus.getPosition();
- List<CalorimeterHit> iHits = iclus.getCalorimeterHits();
- mchits += iHits.size(); // includes hits below threshold
- }
- nclusCheat = mcClusColl.size();
-
- System.out.println(" Comparing cheat clusters: "
- +", cheater="+nclusCheat);
- System.out.println(" Comparing total #hits: "
- +" in event: "+hitmap.size()
- +", in MCClusters="+mchits );
-
- int nhits = hitmap.size();
- AIDA aida = AIDA.defaultInstance();
- aida.cloud1D(colName+"-Nhits").fill( nhits );
- aida.cloud1D(colName+"-NmcHits").fill(mchits);
- aida.cloud1D(colName+"-diffHitsColl-mc").fill(nhits-mchits);
- aida.cloud1D(colName+"-numCheatClusters").fill(nclusCheat);
-
- fillDirectedTreeNtuple( recoClusColl, mcClusColl );
- }
-
- private MCParticle getMCParticleInCluster(BasicCluster clust) {
- for( CalorimeterHit hit : clust.getCalorimeterHits() ) {
- SimCalorimeterHit simhit = (SimCalorimeterHit)hit;
- if( simhit!=null && simhit.getMCParticleCount()==1 ) {
- return simhit.getMCParticle(0);
- }
- }
- return null;
- }
-
- private static Map<String,Readout> _roMap;
- private static RunControlParameters _runPar;
- private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
- private Map<Long,Long> _ParentMap;
- private Map<Long,CalorimeterHit> embhitmap,emehitmap,hdbhitmap,hdehitmap;
- private String _distType,_calType;
- private int _nLyr;
- private int _nZ;
- private int _nPhi;
- private int _evtnum;
-
- private ITree _tree;
- private ITupleFactory _tf;
- private ITuple _tupleEM;
- private ITuple _tupleHD;
- private ITuple _tupleMCem,_tupleMChd;
- private ITuple _tupleEvt;
- private ITuple emClusFolder,emCellFolder;
- private ITuple hdClusFolder,hdCellFolder;
- private ITuple emmcClusFolder,emmcCellFolder;
- private ITuple hdmcClusFolder,hdmcCellFolder;
- private ClusterBuilder clusBuilder = new ClusterBuilder();
-}
CVSspam 0.2.8
