| Commit in lcsim/src/org/lcsim/recon/cluster/density on MAIN | |||
| CalculateDistance.java | +58 | added 1.1 | |
| Cluster.java | +20 | added 1.1 | |
| ClusterBuilder.java | +64 | added 1.1 | |
| DirectedTree.java | +625 | added 1.1 | |
| Kinem.java | +18 | added 1.1 | |
| LoadMyCalorimeterHit.java | +283 | added 1.1 | |
| MCCluster.java | +21 | added 1.1 | |
| MCReconstruction.java | +244 | added 1.1 | |
| MCUtil.java | +33 | added 1.1 | |
| MyCalorimeterHit.java | +34 | added 1.1 | |
| MyTools.java | +89 | added 1.1 | |
| PFlow.java | +113 | added 1.1 | |
| RunControlParameters.java | +191 | added 1.1 | |
| +1793 | |||
GL: first commit for Vishnu's clustering code
diff -N CalculateDistance.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ CalculateDistance.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,58 @@
+package org.lcsim.recon.cluster.density;
+//import java.util.*;
+import Jama.*;
+
+public class CalculateDistance{
+
+ public double CalculateDistance(String distType,double[] a,double[] b){
+
+ double dist = 0.0;
+
+ if(distType=="Euclidean"){
+ double delX = a[0]-b[0];
+ double delY = a[1]-b[1];
+ double delZ = a[2]-b[2];
+ dist = delX*delX+delY*delY+delZ*delZ;
+ if(dist!=0.0) dist = Math.sqrt(dist);
+ }
+
+ if(distType=="Angular"){
+ double theta1 = Kinem.calcTheta(a);
+ double phi1 = Kinem.calcPhi(a);
+ double theta2 = Kinem.calcTheta(b);
+ double phi2 = Kinem.calcPhi(b);
+ double delTheta = Math.abs(theta1-theta2);
+ double delPhi = Math.abs(phi1-phi2);
+ if(delPhi>Math.PI) delPhi = 2.0*Math.PI-delPhi;
+ dist = delTheta*delTheta+delPhi*delPhi;
+ if(dist!=0.0) dist = Math.sqrt(dist);
+ }
+
+ if(distType=="Mahalanobis"){
+ double[][] del = new double[3][1];
+ for(int m=0;m<3;m++){
+ for(int n=0;n<1;n++){
+ del[m][n] = a[m]-b[m];
+ }
+ }
+ Matrix M1 = new Matrix(del);
+ Matrix M2 = M1.transpose();
+
+ double[][] covar = new double[3][3];
+ for(int m=0;m<3;m++){
+ for(int n=0;n<3;n++){
+ covar[m][n] += (a[m]-b[m])*(a[n]-b[n]);
+ }
+ }
+ Matrix cov = new Matrix(covar);
+ Matrix covInv = cov.inverse();
+ Matrix temp1 = M2.times(covInv);
+ Matrix temp2 = temp1.times(M1);
+ dist = temp2.trace();
+ if(dist>0.) dist = Math.sqrt(dist);
+ }
+
+ return dist;
+ }
+
+}
diff -N Cluster.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ Cluster.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,20 @@
+package org.lcsim.recon.cluster.density;
+
+import java.util.Vector;
+
+public class Cluster {
+
+ public Cluster(double ene,double xpos,double ypos,double zpos,
+ double size,Vector<MyCalorimeterHit> hitColl)
+ {
+ _ene = ene;
+ _xpos = xpos;
+ _ypos = ypos;
+ _zpos = zpos;
+ _size = size;
+ _hitColl = hitColl;
+ }
+
+ double _ene,_xpos,_ypos,_zpos,_size;
+ Vector<MyCalorimeterHit> _hitColl;
+}
diff -N ClusterBuilder.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ ClusterBuilder.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,64 @@
+// import hep.aida.*;
+// import hep.lcio.event.*;
+// import hep.lcio.implementation.io.LCFactory;
+// import hep.lcio.io.*;
+// import java.io.IOException;
+// import Jama.Matrix;
+
+package org.lcsim.recon.cluster.density;
+import java.util.*;
+import org.lcsim.event.SimCalorimeterHit;
+
+public class ClusterBuilder {
+
+ public List makeClusters(Vector<Vector<MyCalorimeterHit>> hitColl,RunControlParameters rpar){
+ List<Cluster> clusters = new ArrayList<Cluster>();
+ int nclus = hitColl.size();
+ if(nclus==0) System.out.println("Damn 1");
+ for(int i=0;i<nclus;i++){
+ double cluE = 0.0;
+ double cluX = 0.0; double cluY = 0.0; double cluZ = 0.0;
+ double sumWt = 0.0;
+ double cluSize = hitColl.get(i).size();
+ // System.out.println(cluSize);
+ if(cluSize==0) System.out.println("Damn 2");
+ for( MyCalorimeterHit myhit : hitColl.get(i) ){
+ SimCalorimeterHit hit = myhit.getHit();
+ double ene = hit.getRawEnergy();
+ double[] pos = hit.getPosition();
+ cluE += ene;
+ if(rpar.getCentroidWeightType()=="Energy"){
+ cluX += ene*pos[0];
+ cluY += ene*pos[1];
+ cluZ += ene*pos[2];
+ sumWt += ene;
+ }
+ if(rpar.getCentroidWeightType()=="Density"){
+ if(cluSize>1){
+ double dens = myhit.getDensity(0);
+ cluX += dens*pos[0];
+ cluY += dens*pos[1];
+ cluZ += dens*pos[2];
+ sumWt += dens;
+ }
+ else{
+ cluX = pos[0];
+ cluY = pos[1];
+ cluZ = pos[2];
+ sumWt = 1.0;
+ }
+ }
+ }
+ if(sumWt>0.){
+ cluX = cluX/sumWt;
+ cluY = cluY/sumWt;
+ cluZ = cluZ/sumWt;
+ }
+ // System.out.println(cluE+" "+cluX+" "+cluY+" "+cluZ);
+ Cluster clus = new Cluster(cluE,cluX,cluY,cluZ,cluSize,hitColl.get(i));
+ clusters.add(clus);
+ }
+ return clusters;
+ }
+
+}
diff -N DirectedTree.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ DirectedTree.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,625 @@
+package org.lcsim.recon.cluster.density;
+
+import java.io.IOException;
+import java.util.*;
+
+import hep.aida.ITree;
+import hep.aida.ITupleFactory;
+import hep.aida.ITuple;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimCalorimeterHit;
+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.HitsCluster;
+import hep.physics.vec.Hep3Vector;
+
+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){
+
+ _evtnum = event.getEventNumber();
+ System.out.println("Starting DirectedTree for event no."+" "+_evtnum);
+ _tree.cd("DirectedTree");
+
+ _roMap = event.getDetector().getReadouts();
+ LoadMyCalorimeterHit loader = LoadMyCalorimeterHit.getInstance();
+
+// emhitmap = loader.getEMHitMap(event);
+// hdhitmap = loader.getHDHitMap(event);
+
+ 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";
+// Vector<MyCalorimeterHit> emhit = new Vector<MyCalorimeterHit>();
+// emhit.addAll(embhitmap.values());
+// emhit.addAll(emehitmap.values());
+// emhitmap = new HashMap<Long,MyCalorimeterHit>();
+// emhitmap.putAll(embhitmap);
+// emhitmap.putAll(emehitmap);
+
+ SegmentationBase segm = (SegmentationBase)_roMap.get(embName).getSegmentation();
+ Vector<Vector<MyCalorimeterHit>> embTrees = makeTree(embhitmap, segm);
+ ClusterBuilder embclus = new ClusterBuilder();
+ List<Cluster> embClusters= embclus.makeClusters(embTrees,_runPar);
+
+ List<MCCluster> mcEMB= MCReconstruction.getInstance().getmcClusters(embName);
+
+// List<HitsCluster> embClusters1 = event.get(HitsCluster.class, embName+"Clusters");
+
+ fillDirectedTreeNtuple(embClusters,mcEMB);
+
+
+ _calType = "HD";
+// Vector<MyCalorimeterHit> hdhit = new Vector<MyCalorimeterHit>();
+// hdhit.addAll(hdbhitmap.values());
+// hdhit.addAll(hdehitmap.values());
+// hdhitmap = new HashMap<Long,MyCalorimeterHit>();
+// hdhitmap.putAll(hdbhitmap);
+// hdhitmap.putAll(hdehitmap);
+
+ segm = (SegmentationBase)_roMap.get(hdbName).getSegmentation();
+ Vector<Vector<MyCalorimeterHit>> hdbTrees = makeTree(hdbhitmap, segm);
+ ClusterBuilder hdbclus = new ClusterBuilder();
+ List<Cluster> hdbClusters= hdbclus.makeClusters(hdbTrees,_runPar);
+ List<MCCluster> mcHDB = MCReconstruction.getInstance().getmcClusters(hdbName);
+
+// List<HitsCluster> hdbClusters1 = event.get(HitsCluster.class, hdbName+"Clusters");
+ fillDirectedTreeNtuple(hdbClusters,mcHDB);
+ }
+ 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<MyCalorimeterHit>>
+ makeTree(Map<Long,MyCalorimeterHit> cahitmap, SegmentationBase segm) {
+
+ _ParentMap = new HashMap<Long,Long>();
+ Vector<MyCalorimeterHit> RootVec = new Vector<MyCalorimeterHit>();
+
+ for( MyCalorimeterHit myhit1 : cahitmap.values() ) {
+// for(int i=0;i<cahit.size();i++){
+// MyCalorimeterHit myhit1 = (MyCalorimeterHit) cahit.get(i);
+ SimCalorimeterHit ihit = myhit1.getHit();
+ long cellid = ihit.getCellID();
+ double idens = myhit1.getDensity(0);
+ // 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.){
+ 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];
+ Long jkey = new Long(jid);
+ MyCalorimeterHit myhit2 = cahitmap.get(jkey);
+ if(myhit2==null) continue;
+ SimCalorimeterHit jhit = myhit2.getHit();
+ double jdens = myhit2.getDensity(0);
+ 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(jkey);
+ System.out.println(" jhit: layer="+(jid&0x7f)
+ +", cellID="+MyTools.printID(jid)
+ +", idens="+idens
+ +", jdens="+jdens
+ +", dist="+distance
+ +", densDiff="+densDiff);
+// if(_evtnum==81 && cellid==1672347904) System.out.println(cellid+" "+jid+" "+jkey+" "+myhit2+" "+
+// idens+" "+jdens+" "+densDiff);
+ if(densDiff>maxdensDiff){
+ maxdensDiff = densDiff;
+ maxdensID = jid;
+ }
+ calcD = null;
+ }
+
+// if(_evtnum==81 && cellid== 1672347904) System.out.println("densDiff="+maxdensDiff);
+
+ 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);
+ System.out.println("pt1: parentmap.put: "+MyTools.printID(key1.longValue())+" "+MyTools.printID(key2.longValue()));
+ _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);
+// if(_evtnum==81 && cellid== 1672347904) System.out.println(parent);
+ while(parent!=null){
+ temporary.add(parent);
+ Long lkey = parent;
+ parent = _ParentMap.get(lkey);
+ }
+ if(temporary.contains(key)){
+ removeItems.add(jkey);
+ }
+ }
+ nVec.removeAll(removeItems);
+
+ if(nVec.size()==0){
+ RootVec.add(myhit1);
+ }
+ else{
+ double dmin = 9999.;
+ Long parentKey = new Long(0);
+ for(Long jkey : nVec) {
+// if(_evtnum==81 && cellid== 1672347904) System.out.println(jkey);
+ MyCalorimeterHit myhit2 = 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(key,parentKey);
+ System.out.println("Cell is matched to identical density neighbor");
+ }
+ }
+ }
+ else{
+ RootVec.add(myhit1);
+ System.out.println("Zero density root");
+ }
+ }
+
+ Vector<MyCalorimeterHit> startingPoints = new Vector<MyCalorimeterHit>();
+ for( MyCalorimeterHit myhit1 : cahitmap.values() ) {
+// for(int i=0;i<cahit.size();i++){
+// MyCalorimeterHit myhit1 = (MyCalorimeterHit) cahit.get(i);
+ SimCalorimeterHit ihit = myhit1.getHit();
+ long cellid = ihit.getCellID();
+ Long value = new Long(cellid);
+ boolean isParent = _ParentMap.containsValue(value);
+ boolean isRoot = RootVec.contains(myhit1);
+ if(!isParent && !isRoot){
+ startingPoints.add(myhit1);
+ }
+ }
+
+ Vector<Vector<MyCalorimeterHit>> branches
+ = new Vector<Vector<MyCalorimeterHit>>();
+ for(int i=0;i<startingPoints.size();i++) {
+ branches.add( new Vector<MyCalorimeterHit>() );
+ }
+ for(int i=0;i<startingPoints.size();i++){
+ MyCalorimeterHit myhit1 = startingPoints.get(i);
+ branches.get(i).add(myhit1);
+ SimCalorimeterHit ihit = myhit1.getHit();
+ long cellid = ihit.getCellID();
+ Long key = new Long(cellid);
+ while(_ParentMap.containsKey(key)){
+ Long jkey = _ParentMap.get(key);
+ MyCalorimeterHit myhit2 = cahitmap.get(jkey);
+ branches.get(i).add(myhit2);
+// if(_evtnum==81) System.out.println("filling branches"+" "+i+" "+myhit1+" "+key+" "+jkey+" "+myhit2);
+ key = jkey;
+ }
+ }
+
+ int vsiz = RootVec.size();
+ // System.out.println("no. of roots"+" "+vsiz);
+ Vector<Vector<MyCalorimeterHit>> caTrees = new Vector<Vector<MyCalorimeterHit>>();
+ for(int i=0;i<vsiz;i++){
+ caTrees.add( new Vector<MyCalorimeterHit>() );
+ }
+
+ for(int i=0;i<vsiz;i++){
+ MyCalorimeterHit myhit1 = RootVec.get(i);
+ long tmpID = myhit1.getHit().getCellID();
+ System.out.println("ROOT id="+MyTools.printID(tmpID));
+ caTrees.get(i).add(myhit1);
+ for(int j=0;j<startingPoints.size();j++) {
+ MyCalorimeterHit parentHit = branches.get(j).lastElement();
+// if(_evtnum==81) System.out.println(myhit1+" "+parentHit+" "+startingPoints.size()+" "+j);
+ if(parentHit.equals(myhit1)) {
+ for( MyCalorimeterHit myhit2 : branches.get(j) ) {
+ if(!caTrees.get(i).contains(myhit2)){
+ caTrees.get(i).add(myhit2);
+ long tmppID = myhit2.getHit().getCellID();
+ System.out.println("Branches id="+MyTools.printID(tmppID));
+ }
+ }
+ }
+ }
+ }
+
+ 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<MCCluster> mctruth)
+ {
+ System.out.println("fillDirectedTreeNtuple() called, "+recon+" "+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<MyCalorimeterHit,Integer> recTagMap
+ = new HashMap<MyCalorimeterHit,Integer>();
+ for(Cluster clust : recon ) {
+ Vector<MyCalorimeterHit> cellVec = clust._hitColl;
+ // 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(MyCalorimeterHit myhit : cellVec) {
+ SimCalorimeterHit ihit = myhit.getHit();
+ Integer tagObject = new Integer(tag);
+ System.out.println("putting <"+myhit+"->"+tagObject+">");
+ recTagMap.put(myhit,tagObject);
+ double[] pos = ihit.getPosition();
+ double cdens = myhit.getDensity(0);
+ 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,(int) clust._size);
+ 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,(int) clust._size);
+ hdClusFolder.addRow();
+ }
+ tag++;
+ }
+
+ tag = 0;
+ for( MCCluster clust : mctruth ) {
+ Vector<MyCalorimeterHit> cellVec = clust._hitColl;
+ for(MyCalorimeterHit myhit : cellVec) {
+ Integer tagObject = recTagMap.get(myhit);
+ System.out.println("getting: "+myhit+" -> "+tagObject);
+ int rtag = -1;
+ if(tagObject!=null) rtag = tagObject.intValue();
+ SimCalorimeterHit ihit = myhit.getHit();
+ double[] pos = ihit.getPosition();
+ double cdens = myhit.getDensity(0);
+ 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();
+ // if(tag==17){
+ System.out.println("EM Clus cells tag="+tag
+ +", 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 = clust._mcpart;
+ float charge = (float)mcpart.getCharge();
+ System.out.println("charge="+charge);
+ int particleID = mcpart.getType().getPDGID();
+ float mcE = (float)mcpart.getEnergy();
+ Hep3Vector mcMom = mcpart.getMomentum();
+ if(_calType=="EM"){
+ emmcClusFolder.fill(0,(float) clust._ene);
+ emmcClusFolder.fill(1,(float) clust._xpos);
+ emmcClusFolder.fill(2,(float) clust._ypos);
+ emmcClusFolder.fill(3,(float) clust._zpos);
+ emmcClusFolder.fill(4,(int) clust._size);
+ 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._ene);
+ hdmcClusFolder.fill(1,(float) clust._xpos);
+ hdmcClusFolder.fill(2,(float) clust._ypos);
+ hdmcClusFolder.fill(3,(float) clust._zpos);
+ hdmcClusFolder.fill(4,(int) clust._size);
+ 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();
+
+ }
+
+// public SegmentationBase getSegmentation(long cellid) {
+// // warning: this is very detector specific!!!
+// int segmID = cellid & 0x0;
+// }
+
+ private static Map<String,Readout> _roMap;
+ private static RunControlParameters _runPar;
+ private Map<Long,Long> _ParentMap;
+ private Map<Long,MyCalorimeterHit> emhitmap,hdhitmap,cahitmap;
+ private Map<Long,MyCalorimeterHit> 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;
+}
diff -N Kinem.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ Kinem.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,18 @@
+package org.lcsim.recon.cluster.density;
+
+public class Kinem {
+
+ public static double calcTheta(double[] pos){
+ double theta = 0.0;
+ double r = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
+ theta = Math.atan2(r,pos[2]);
+ return theta;
+ }
+
+ public static double calcPhi(double[] pos){
+ double phi = 0.0;
+ phi = Math.atan2(pos[1],pos[0]);
+ if(phi<0.0) phi += Math.PI;
+ return phi;
+ }
+}
diff -N LoadMyCalorimeterHit.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ LoadMyCalorimeterHit.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,283 @@
+package org.lcsim.recon.cluster.density;
+
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Map;
+import java.util.HashMap;
+
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.geometry.segmentation.SegmentationBase;
+import org.lcsim.geometry.compact.Readout;
+
+public class LoadMyCalorimeterHit {
+
+ public static LoadMyCalorimeterHit getInstance(){
+ if(_me==null){
+ _me = new LoadMyCalorimeterHit();
+ }
+ return _me;
+ }
+
+ public void process(EventHeader event){
+ System.out.println("Running LoadMyCalorimeterHit.process");
+ this.reset();
+ _segmMap = event.getDetector().getReadouts();
+ _runPar = RunControlParameters.getInstance();
+ _event = event;
+ }
+
+ private LoadMyCalorimeterHit(){
+ if(_me==null) _me=this;
+ _collMap = new HashMap<String,Map<Long,MyCalorimeterHit>>();
+ _myEMhits = new HashMap<Long,MyCalorimeterHit>();
+ _myHDhits = new HashMap<Long,MyCalorimeterHit>();
+ }
+
+ public void reset(){
+ System.out.println("reset: EM="+_myEMhits.size()+", HD="+_myHDhits.size());
+ _myEMhits.clear();
+ _myHDhits.clear();
+ System.out.println("reset: EM="+_myEMhits.size()+", HD="+_myHDhits.size());
+ // Loop over all existing hit collections
+ for( String key : _collMap.keySet() ) {
+ // clear hit maps for each subdetector collection
+ Map<Long,MyCalorimeterHit> subDetHitMap = _collMap.get(key);
+ System.out.println("reset: key="+key+", size="+subDetHitMap.size());
+ subDetHitMap.clear();
+ System.out.println("reset: key="+key+", size="+subDetHitMap.size());
+ }
+ System.out.println("reset: #colls="+_collMap.size());
+ _collMap.clear();
+ System.out.println("reset: #colls="+_collMap.size());
+ }
+
+// public Map<Long,MyCalorimeterHit> getEMHitMap(EventHeader event){
+// if(_myEMhits.size()==0){
+// List<SimCalorimeterHit> ecol = new ArrayList<SimCalorimeterHit>();
+// ecol.addAll( event.getSimCalorimeterHits("EcalBarrHits") );
+// ecol.addAll( event.getSimCalorimeterHits("EcalEndcapHits") );
+// System.out.println("getEMHitMap: sizes="
+// +event.getSimCalorimeterHits("EcalBarrHits").size()
+// +" "+event.getSimCalorimeterHits("EcalEndcapHits").size()
+// +", total="+ecol.size());
+// // System.out.println(ecol);
+// String calType = "EM";
+// fillHitMap(ecol,_myEMhits,calType);
+// }
+// return _myEMhits;
+// }
+
+// public Map<Long,MyCalorimeterHit> getHDHitMap(EventHeader event){
+// if(_myHDhits.size()==0){
+// List<SimCalorimeterHit> hcol = new ArrayList<SimCalorimeterHit>();
+// hcol.addAll( event.getSimCalorimeterHits("HcalBarrHits") );
+// hcol.addAll( event.getSimCalorimeterHits("HcalEndcapHits") );
+// System.out.println("getHDHitMap: sizes="
+// +event.getSimCalorimeterHits("HcalBarrHits").size()
+// +" "+event.getSimCalorimeterHits("HcalEndcapHits").size()
+// +", total="+hcol.size());
+
+// String calType = "HD";
+// fillHitMap(hcol,_myHDhits,calType);
+// }
+// return _myHDhits;
+// }
+
+ /**
+ * Returns a CellID decoder
+ * @param colName name of collection with data to be returned
+ */
+ public final Map<Long,MyCalorimeterHit> getCollHitMap(final String colName) {
+ Map<Long,MyCalorimeterHit> retColl = _collMap.get( colName );
+ if(retColl==null) {
+ // if it does not exist, create it now
+ retColl = new HashMap<Long,MyCalorimeterHit>();
+ _collMap.put(colName, retColl);
+ }
+ assert retColl!=null : "retColl is empty for colName="+colName;
+ String calType = "??";
+ if( colName.contains("Ecal") ) calType = "EM";
+ if( colName.contains("Hcal") ) calType = "HD";
+ assert !calType.contains("??") : "Unknown calorimeter type.";
+
+ System.out.println("getCollHitMap: name=<"+colName
+ +">, size="+retColl.size()
+ +", calType="+calType);
+ if( retColl.size()==0 ) fillHitMap( colName, retColl, calType );
+ System.out.println("getCollHitMap: name=<"+colName
+ +">, size="+retColl.size());
+ return retColl;
+ }
+
+ /**
+ * Returns a hit map, with all hits contained in a given collection.
+ * @param colName name of collection with data to be returned
+ */
+ public final SegmentationBase getSegmentation(final String colName) {
+ return (SegmentationBase)_segmMap.get(colName).getSegmentation();
+ }
+
+ public void fillHitMap(final String colName,
+ Map<Long,MyCalorimeterHit> hitmap,
+ String calType)
+ {
+ if(calType=="EM"){
+ _weight = _runPar.getEMweight();
+ _mip = _runPar.getEMmip();
+ _thresh = _runPar.getEMthresh();
+ _nLyr = _runPar.getLyrNeighEM();
+ _nZ = _runPar.getZNeighEM();
+ _nPhi = _runPar.getPhiNeighEM();
+ }
+ else {
+ _weight = _runPar.getHDweight();
+ _mip = _runPar.getHDmip();
+ _thresh = _runPar.getHDthresh();
+ _nLyr = _runPar.getLyrNeighHD();
+ _nZ = _runPar.getZNeighHD();
+ _nPhi = _runPar.getPhiNeighHD();
+ }
+ /*
+ _nLyr = _runPar.getLyrNeigh1();
+ _nZ = _runPar.getZNeigh1();
+ _nPhi = _runPar.getPhiNeigh1();
+ */
+
+ try {
+ List<SimCalorimeterHit> tmpHits= _event.getSimCalorimeterHits(colName);
+ for( SimCalorimeterHit ihit : tmpHits ){
+ double ene = _weight*(ihit.getRawEnergy());
+ if(ene>_thresh*_mip){
+ long cellid = ihit.getCellID();
+ Long key = new Long(cellid);
+ MyCalorimeterHit khit = new MyCalorimeterHit(ihit);
+ hitmap.put(key,khit);
+ }
+ }
+
+ for( MyCalorimeterHit khit : hitmap.values() ) {
+ SimCalorimeterHit ihit = khit.getHit();
+ long cellid = ihit.getCellID();
+ // System.out.println("cell index (lyr,theta,phi)"+" "+MyTools.getCalLayer(cellid)+
+ // " "+MyTools.getThetaBin(cellid)+" "+MyTools.getPhiBin(cellid));
+
+ // int neigh[] = _cPar.getNeighbors(cellid,_nLyr,_nZ,_nPhi);
+
+ SegmentationBase segm = (SegmentationBase)_segmMap.get(colName).getSegmentation();
+ segm.setID(cellid);
+
+ long neigh[] = segm.getNeighbourIDs(_nLyr,_nZ,_nPhi);
+ System.out.println(" Density: id="+MyTools.printID(cellid)
+ +", #neighbors="+neigh.length);
+ double nNeigh = 0.0;
+ for(int j=0;j<neigh.length;j++){
+ long jid = neigh[j];
+ // System.out.println("neigh index"+" "+MyTools.getCalLayer(jid)+" "+
+ // MyTools.getThetaBin(jid)+" "+MyTools.getPhiBin(jid));
+ Long jkey = new Long(jid);
+ MyCalorimeterHit lhit = (MyCalorimeterHit) hitmap.get(jkey);
+ if(lhit==null) continue;
+ // System.out.println("not null");
+// jhit = lhit.getHit();
+ nNeigh++;
+ }
+ khit.setDensity(0,nNeigh);
+ khit.setDensity(1,0.0);
+ AIDA.defaultInstance().cloud1D("density").fill(nNeigh);
+ }
+ }
+ catch(Exception e) {
+ System.out.println(" Evt "+_event.getEventNumber()
+ +": No data in collection "+colName);
+ // This happens quite frequently, no data from this component
+ }
+ }
+
+ public void fillHitMap(final List<SimCalorimeterHit> col,
+ Map<Long,MyCalorimeterHit> hitmap,
+ String calType)
+ {
+ if(calType=="EM"){
+ _weight = _runPar.getEMweight();
+ _mip = _runPar.getEMmip();
+ _thresh = _runPar.getEMthresh();
+ _nLyr = _runPar.getLyrNeighEM();
+ _nZ = _runPar.getZNeighEM();
+ _nPhi = _runPar.getPhiNeighEM();
+ }
+ else {
+ _weight = _runPar.getHDweight();
+ _mip = _runPar.getHDmip();
+ _thresh = _runPar.getHDthresh();
+ _nLyr = _runPar.getLyrNeighHD();
+ _nZ = _runPar.getZNeighHD();
+ _nPhi = _runPar.getPhiNeighHD();
+ }
+ /*
+ _nLyr = _runPar.getLyrNeigh1();
+ _nZ = _runPar.getZNeigh1();
+ _nPhi = _runPar.getPhiNeigh1();
+ */
+
+ // System.out.println(col.getNumberOfElements());
+// for(int i=0;i<col.size();i++){
+// ihit = (SimCalorimeterHit) col.getElementAt(i);
+ for( SimCalorimeterHit ihit : col ){
+ double ene = _weight*(ihit.getRawEnergy());
+ if(ene>_thresh*_mip){
+ long cellid = ihit.getCellID();
+ Long key = new Long(cellid);
+ MyCalorimeterHit khit = new MyCalorimeterHit(ihit);
+ hitmap.put(key,khit);
+ }
+ }
+
+ for( MyCalorimeterHit khit : hitmap.values() ) {
+ SimCalorimeterHit ihit = khit.getHit();
+ long cellid = ihit.getCellID();
+ // System.out.println("cell index (lyr,theta,phi)"+" "+MyTools.getCalLayer(cellid)+
+ // " "+MyTools.getThetaBin(cellid)+" "+MyTools.getPhiBin(cellid));
+
+ // int neigh[] = _cPar.getNeighbors(cellid,_nLyr,_nZ,_nPhi);
+
+ SegmentationBase segm = (SegmentationBase)_segmMap.get("EcalBarrHits").getSegmentation();
+ segm.setID(cellid);
+
+ long neigh[] = segm.getNeighbourIDs(_nLyr,_nZ,_nPhi);
+ double nNeigh = 0.0;
+ for(int j=0;j<neigh.length;j++){
+ long jid = neigh[j];
+ // System.out.println("neigh index"+" "+MyTools.getCalLayer(jid)+" "+
+ // MyTools.getThetaBin(jid)+" "+MyTools.getPhiBin(jid));
+ Long jkey = new Long(jid);
+ MyCalorimeterHit lhit = (MyCalorimeterHit) hitmap.get(jkey);
+ if(lhit==null) continue;
+ // System.out.println("not null");
+// jhit = lhit.getHit();
+ nNeigh++;
+ }
+ khit.setDensity(0,nNeigh);
+ khit.setDensity(1,0.0);
+ AIDA.defaultInstance().cloud1D("density").fill(nNeigh);
+ }
+ }
+
+
+ private static LoadMyCalorimeterHit _me = null;
+ /** Current event */
+ private EventHeader _event;
+ private static RunControlParameters _runPar;
+ private double _weight;
+ private double _mip;
+ private double _thresh;
+ private int _nLyr;
+ private int _nZ;
+ private int _nPhi;
+ private Map<String,Map<Long,MyCalorimeterHit>> _collMap;
+ private Map<Long,MyCalorimeterHit> _myEMhits;
+ private Map<Long,MyCalorimeterHit> _myHDhits;
+ private Map<String,Readout> _segmMap;
+}
diff -N MCCluster.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MCCluster.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,21 @@
+package org.lcsim.recon.cluster.density;
+//import hep.lcio.event.*;
+import java.util.Vector;
+import org.lcsim.event.MCParticle;
+
+public class MCCluster {
+
+ public MCCluster(MCParticle mcpart,double ene,double xpos,double ypos,double zpos,double size,Vector<MyCalorimeterHit> hitColl){
+ _mcpart = mcpart;
+ _ene = ene;
+ _xpos = xpos;
+ _ypos = ypos;
+ _zpos = zpos;
+ _size = size;
+ _hitColl = hitColl;
+ }
+
+ double _ene,_xpos,_ypos,_zpos,_size;
+ MCParticle _mcpart;
+ Vector<MyCalorimeterHit> _hitColl;
+}
diff -N MCReconstruction.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MCReconstruction.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,244 @@
+package org.lcsim.recon.cluster.density;
+
+import java.io.IOException;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+
+public class MCReconstruction extends Driver {
+
+ static MCReconstruction getInstance(){
+ if(_me==null){
+ _me = new MCReconstruction();
+ }
+ return _me;
+ }
+
+ public void process(EventHeader event){
+ System.out.println("Running _mcreco.process");
+ this.reset();
+// _cPar = CalParameters.getInstance(event);
+ _runPar = RunControlParameters.getInstance();
+ _event = event;
+ }
+
+ private MCReconstruction(){
+ if(_me==null) _me=this;
+ _myEMmcClusters = new ArrayList<MCCluster>();
+ _myHDmcClusters = new ArrayList<MCCluster>();
+ }
+
+ public void reset(){
+ _myEMmcClusters.clear();
+ _myHDmcClusters.clear();
+ _mcClusters.clear();
+ }
+
+ public List<MCCluster> getmcClusters(final String colName) {
+ if(_mcClusters.size()==0) {
+ Map<Long,MyCalorimeterHit> hitmap
+ = LoadMyCalorimeterHit.getInstance().getCollHitMap(colName);
+ List<MCParticle> mcp = _event.getMCParticles();
+ fillMCclusterList(hitmap, mcp, _mcClusters);
+ }
+ return _mcClusters;
+ }
+
+ public void fillMCclusterList(Map<Long,MyCalorimeterHit> hitmap,
+ List<MCParticle> mcpcol,
+ List<MCCluster> mclus)
+ {
+
+ Vector<MCParticle> mcpVec = new Vector<MCParticle>();
+// for(int i=0;i<mcpcol.size();++i){
+// MCParticle mcp = (MCParticle) mcpcol.getElementAt(i);
+ for( MCParticle mcp : mcpcol ) {
+/// This portion needs to be clarified with Guilherme
+// boolean simPart = mcp.isCreatedInSimulation();
+ // if(!simPart){
+ // mcpVec.add(mcp);
+ // }
+
+ boolean preCalorimeter = MCUtil.createdBeforeCal(mcp);
+ boolean calWorthy = MCUtil.isCalorimeterRelevant(mcp);
+ int genStatus = mcp.getGeneratorStatus();
+ if(preCalorimeter && calWorthy && genStatus<3){
+ mcpVec.add(mcp);
+ }
+ }
+
+ int nmcp = mcpVec.size();
+ Map<MCParticle, Vector<MyCalorimeterHit>> mcpHit
+ = new HashMap<MCParticle, Vector<MyCalorimeterHit>>();
+ for( MCParticle mcp : mcpVec ) {
+ System.out.println("Add hits vector for "+mcp);
+ mcpHit.put( mcp, new Vector<MyCalorimeterHit>() );
+ }
+
+ for( MyCalorimeterHit myhit : hitmap.values() ) {
+ SimCalorimeterHit ihit = myhit.getHit();
+ int nmc = ihit.getMCParticleCount();
+ MCParticle parent = null;
+ for(int imc=0;imc<nmc;imc++){
+ MCParticle ipart = ihit.getMCParticle(imc);
+ int ndx = mcpVec.indexOf(ipart);
+ if(ndx<0){
+ int nparent = ipart.getParents().size();
+ int j = 0; boolean traceBack = true;
+ while(traceBack && j<nparent){
+ parent = ipart.getParents().get(j);
+ boolean preCalorimeter = MCUtil.createdBeforeCal(parent);
+
+// System.out.println("starting traceback"+" "+ipart.getPDG()
+// +" "+ipart.getGeneratorStatus()
+// +" "+ipart.isCreatedInSimulation()
+// +" "+MCUtil.createdBeforeCal(ipart)
+// +" "+MCUtil.isCalorimeterRelevant(ipart)
+// +" "+nparent+" "+j);
+
+ if(preCalorimeter){
+ ndx = mcpVec.indexOf(parent);
+ ipart = parent;
+// System.out.println("in traceback"+" "+ndx+" "+parent.getPDG()
+// +" "+parent.getGeneratorStatus());
+ traceBack = false;
+ }
+ j++;
+ }
+ }
+ if(ndx>=0) {
+ mcpHit.get(ipart).add(myhit);
+ }
+ }
+ }
+
+ makeMCclusters(mcpVec,mcpHit,mclus);
+ }
+
+ public void makeMCclusters(Vector<MCParticle> mcp,
+ Map<MCParticle, Vector<MyCalorimeterHit>> hitColl,
+ List<MCCluster> mcclusters)
+ {
+ int nclus = hitColl.size();
+ for( MCParticle mcpart : mcp ) {
+ double cluE = 0.0;
+ double cluX = 0.0; double cluY = 0.0; double cluZ = 0.0;
+ double sumWt = 0.0;
+ Vector<MyCalorimeterHit> cluHits = hitColl.get(mcpart);
+ double cluSize = cluHits.size();
+
+ // if(cluSize==0){
+ /*
+ float[] mom = mcpart.getMomentum();
+ double[] ep = mcpart.getEndpoint();
+ System.out.println("mcpart"+" "+mcpart.getPDG()+" "+mcpart.getGeneratorStatus()+" "+
+ mcpart.isDecayedInCalorimeter()+" "+mom[0]+" "+mom[1]+
+ " "+mom[2]+" "+ep[0]+" "+ep[1]+" "+ep[2]+" "+
+ mcpart.isStopped()+" "+cluSize+" "+i);
+ */
+ // }
+
+ if(cluSize==0) continue;
+ Iterator itr = cluHits.iterator();
+ while(itr.hasNext()){
+ MyCalorimeterHit myhit = (MyCalorimeterHit)itr.next();
+ SimCalorimeterHit hit = myhit.getHit();
+ int npart = hit.getMCParticleCount();
+ double ene = 0;
+ for(int j=0;j<npart;j++){
+ MCParticle jpart = hit.getMCParticle(j);
+ if(jpart.equals(mcpart)){
+ ene = hit.getContributedEnergy(j);
+ }
+ }
+ double[] pos = hit.getPosition();
+ cluE += (double) ene;
+ if(_runPar.getCentroidWeightType()=="Energy"){
+ cluX += (double) ene*pos[0];
+ cluY += (double) ene*pos[1];
+ cluZ += (double) ene*pos[2];
+ sumWt += (double) ene;
+ }
+ if(_runPar.getCentroidWeightType()=="Density"){
+ double dens = myhit.getDensity(0);
+//density sharing needs to be worked out
+ cluX += (double) dens*pos[0];
+ cluY += (double) dens*pos[1];
+ cluZ += (double) dens*pos[2];
+ sumWt += dens;
+ }
+ }
+ if(sumWt>0.){
+ cluX = cluX/sumWt;
+ cluY = cluY/sumWt;
+ cluZ = cluZ/sumWt;
+ }
+ // System.out.println(cluE);
+
+ if(_runPar.Prune()){
+ double cluXorig = cluX;
+ double cluYorig = cluY;
+ double cluZorig = cluZ;
+
+ cluE = 0.0; sumWt = 0.0;
+ cluX = 0.0; cluY = 0.0; cluZ = 0.0;
+ Iterator itr1 = cluHits.iterator();
+ while(itr1.hasNext()){
+ MyCalorimeterHit myhit = (MyCalorimeterHit) itr1.next();
+ SimCalorimeterHit hit = myhit.getHit();
+ int npart = hit.getMCParticleCount();
+ double ene = 0;
+ for(int j=0;j<npart;j++){
+ MCParticle jpart = hit.getMCParticle(j);
+ if(jpart.equals(mcpart)){
+ ene = hit.getContributedEnergy(j);
+ }
+ }
+ double[] pos = hit.getPosition();
+ double diffX = pos[0]-cluXorig;
+ double diffY = pos[1]-cluYorig;
+ double diffZ = pos[2]-cluZorig;
+ double distance = diffX*diffX+diffY*diffY+diffZ*diffZ;
+ if(distance>0.0) distance = Math.sqrt(distance);
+ if(distance<_runPar.getPruningDist()){
+ cluE += (double) ene;
+ if(_runPar.getCentroidWeightType()=="Energy"){
+ cluX += (double) ene*pos[0];
+ cluY += (double) ene*pos[1];
+ cluZ += (double) ene*pos[2];
+ sumWt += (double) ene;
+ }
+ if(_runPar.getCentroidWeightType()=="Density"){
+ double dens = myhit.getDensity(0);
+//density sharing needs to be worked out
+ cluX += (double) dens*pos[0];
+ cluY += (double) dens*pos[1];
+ cluZ += (double) dens*pos[2];
+ sumWt += dens;
+ }
+ }
+ }
+ if(sumWt>0.){
+ cluX = cluX/sumWt;
+ cluY = cluY/sumWt;
+ cluZ = cluZ/sumWt;
+ }
+ }
+
+ // System.out.println(cluE+" "+cluSize);
+ MCCluster clus = new MCCluster(mcpart,cluE,cluX,cluY,cluZ,cluSize,cluHits);
+ mcclusters.add(clus);
+ }
+ }
+
+ private static MCReconstruction _me = null;
+// private static CalParameters _cPar;
+ private EventHeader _event;
+ private static RunControlParameters _runPar;
+ private Map<Long,MyCalorimeterHit> emhitmap,hdhitmap;
+ private static List<MCCluster> _myEMmcClusters = new ArrayList<MCCluster>();
+ private static List<MCCluster> _myHDmcClusters = new ArrayList<MCCluster>();
+ private static List<MCCluster> _mcClusters = new ArrayList<MCCluster>();
+}
diff -N MCUtil.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MCUtil.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,33 @@
+package org.lcsim.recon.cluster.density;
+
+import org.lcsim.event.MCParticle;
+import hep.physics.vec.Hep3Vector;
+
+public class MCUtil{
+
+ public static boolean createdBeforeCal(MCParticle mcpart){
+ boolean beforeCalorimeter = false;
+ Hep3Vector vertex = mcpart.getOrigin();
+ double[] vert = new double[3];
+ vert[0] = vertex.x();
+ vert[1] = vertex.y();
+ vert[2] = vertex.z();
+// for(int ii=0;ii<3;ii++){
+// vert = vertex/10.;
+// }
+ double rad = vert[0]*vert[0]+vert[1]*vert[1];
+ if(rad>0.) rad = Math.sqrt(rad);
+ if(rad<125. && Math.abs(vert[2])<185.){
+ beforeCalorimeter = true;
+ }
+ return beforeCalorimeter;
+ }
+
+ public static boolean isCalorimeterRelevant(MCParticle mcpart){
+ boolean interesting = true;
+ int pid = Math.abs(mcpart.getType().getPDGID());
+ if(pid==12 || pid==14 || pid==16) interesting = false;
+ // if(mcpart.getEnergy()<0.5) interesting = false;
+ return interesting;
+ }
+}
diff -N MyCalorimeterHit.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyCalorimeterHit.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,34 @@
+package org.lcsim.recon.cluster.density;
+
+import org.lcsim.event.SimCalorimeterHit;
+
+public class MyCalorimeterHit {
+
+ public MyCalorimeterHit(SimCalorimeterHit myHit){
+ _myHit = myHit;
+ initDensity();
+ }
+
+ private void initDensity(){
+ _dens = new double[2];
+ _dens[0] = 0.0;
+ _dens[1] = 0.0;
+ }
+
+ public void setDensity(int index, double density){
+ if(index>=_dens.length) throw new IllegalArgumentException();
+ _dens[index] = density;
+ }
+
+ public double getDensity(int index){
+ if(index>=_dens.length) throw new IllegalArgumentException();
+ return this._dens[index];
+ }
+
+ public SimCalorimeterHit getHit(){
+ return this._myHit;
+ }
+
+ private SimCalorimeterHit _myHit;
+ private double [] _dens;
+}
diff -N MyTools.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyTools.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,89 @@
+/*
+ * MyTools.java
+ *
+ * Created on Nov. 2003
+ * G.Lima, J. McCormick (NICADD/NIU)
+ */
+
+package org.lcsim.recon.cluster.density;
+import org.lcsim.event.EventHeader;
+
+/**
+ *
+ * @author G.Lima, J.McCormick
+ */
+public class MyTools
+{
+// // This returns the index of object obj in collection col
+// public static int indexOf( LCObject obj, LCCollection col )
+// {
+// int fail = -1;
+// if( obj == null ) return fail;
+
+// for(int i=0; i<col.getNumberOfElements(); i++) {
+// if( obj == col.getElementAt(i) ) return i;
+// }
+
+// return fail;
+// }
+
+// public static void listCollections(EventHeader evt)
+// {
+// String[] colls = evt.getCollectionNames();
+// for(int i=0; i<colls.length; i++) {
+// System.out.println("Collection name: "+colls[i]);
+// }
+// }
+
+ public static int getSysno( long cellid ) {
+ int sysno = (int)(cellid>>16)&0x7;
+ return sysno;
+ }
+
+ public static int getSysBE( long cellid ) {
+ return (int)((cellid>>16) & 0x3f);
+ }
+
+ public static int getLayer( long cellid ) {
+ int layer = (int)((cellid>>0) & 0xffff);
+ return layer;
+ }
+
+ public static int getThetaBin( long cellid ) {
+ int max = 1<<16; // for 11 bits
+ int thetabin = (int)((cellid>>48) & 0xffff);
+ if( thetabin > max/2-1 ) thetabin = thetabin - max;
+ return thetabin;
+ }
+
+ public static int getPhiBin( long cellid ) {
+// int max = 1<<11; // for 11 bits
+// int phibin = (int)((cellid>>43) & 0x7ff);
+// if( phibin > max/2-1 ) phibin = phibin - max;
+// return phibin;
+ return (int)((cellid>>32) & 0xffff);
+ }
+
+ public static int getECalLayerTB( float[] pos ) {
+ double aux = (584.30 - pos[1])/5;
+ int layer = (int)(aux+0.5);
+ if(layer>29 || layer<0) System.out.println(" Bad ECal layer="+layer);
+ return layer;
+ }
+
+ public static int getHCalLayerTB( float[] pos ) {
+ double aux = (435.0 - pos[1])/25;
+ int layer = (int)(aux+0.5);
+ if(layer>34 || layer<0) System.out.println(" Bad HCal layer="+layer+", pos="+pos[1]);
+ return (int)(aux+0.5);
+ }
+
+ public static String printID(long id) {
+ return new String("<"+Long.toHexString(id)
+ +": "+Integer.toString( getSysBE(id) )
+ +" "+Integer.toString( getLayer(id) )
+ +" "+Integer.toString( getThetaBin(id))
+ +" "+Integer.toString( getPhiBin(id) )
+ +">");
+ }
+}
diff -N PFlow.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ PFlow.java 12 Aug 2005 19:59:49 -0000 1.1 @@ -0,0 +1,113 @@
+
+package org.lcsim.recon.cluster.density;
+
+import hep.aida.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.EventHeader;
+import java.io.File;
+import java.io.IOException;
+
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+import org.lcsim.util.loop.LCSimLoop;
+import org.lcsim.recon.cluster.cheat.CheatClusterDriver;
+
+public class PFlow extends Driver {
+
+ public static void main(String [] args) throws Exception{
+
+ if(args.length==0) help();
+
+ // Setup LCSim event loop
+ LCSimLoop loop = new LCSimLoop();
+
+ // setup input/outfile files
+ File input = new File(args[0]);
+ loop.setLCIORecordSource(input);
+
+// ILCFactory factory = LCFactory.getInstance();
+// LCReader reader = factory.createLCReader();
+// reader.open(args[0]);
+
+ // drivers needed
+ PFlow pflow = new PFlow();
+ loop.add( pflow );
+
+// // ModeDetection anal2 = new ModeDetection(_tf,_tree);
+// // Canny anal3 = new Canny(_tf,_tree);
+
+// reader.registerLCEventListener(mycalhit);
+// reader.registerLCEventListener(mcreco);
+// reader.registerLCEventListener(anal1);
+// // reader.registerLCEventListener(anal2);
+// // reader.registerLCEventListener(anal3);
+
+ // last driver, if output is needed
+// loop.add( new LCIODriver( new File("output.slcio") );
+
+ System.out.println("Starting loop over events");
+ loop.loop(10); // -1 for all
+
+ loop.dispose();
+ _tree.commit();
+
+ org.lcsim.util.aida.AIDA.defaultInstance().saveAs("myhistos.aida");
+ System.out.println("Analyzed "+ pflow._nevt +" events.");
+ }
+
+ public void process(final EventHeader evt) {
+ // event counter
+ ++_nevt;
+
+ // subdrivers processing
+ System.out.print("### Event "+_nevt+" ");
+
+ System.out.print("..calhit..");
+ _calhit.process(evt);
+// System.out.print("..cheater..");
+// _cheater.process(evt);
+
+ System.out.print("..mcreco..");
+ MCReconstruction.getInstance().process(evt);
+ System.out.print("..anal..");
+ _anal.process(evt);
+ System.out.println("..Done.");
+ }
+
+ private static void help(){
+ System.out.println("Usage:");
+ System.out.println(" java "+PFlow.class.getName()+" <input-file>");
+ System.exit(1);
+ }
+
+ public PFlow() throws IOException{
+ _calhit = LoadMyCalorimeterHit.getInstance();
+// _cheater = new CheatClusterDriver();
+
+ prepareNtupleFactory();
+ _anal = new DirectedTree(_tf,_tree);
+
+// add(_cheater);
+ add( MCReconstruction.getInstance() );
+ add(_anal);
+ }
+
+ public void prepareNtupleFactory() throws IOException{
+ _aida = AIDA.defaultInstance();
+
+// _tree=af.createTreeFactory().create("canny-test2.aida","xml",false,true);
+ _tree = _aida.tree();
+ _tf = _aida.analysisFactory().createTupleFactory(_tree);
+ }
+
+ //=== FIELDS ====
+
+ private static AIDA _aida;
+ private static ITree _tree;
+ private static ITupleFactory _tf;
+
+ private int _nevt = 0;
+ private LoadMyCalorimeterHit _calhit;
+ private CheatClusterDriver _cheater;
+ private DirectedTree _anal;
+}
diff -N RunControlParameters.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ RunControlParameters.java 12 Aug 2005 19:59:50 -0000 1.1 @@ -0,0 +1,191 @@
+package org.lcsim.recon.cluster.density;
+
+public class RunControlParameters {
+
+ static RunControlParameters getInstance() {
+ if(_me==null) {
+ _me = new RunControlParameters();
+ }
+ return _me;
+ }
+
+ private RunControlParameters() {
+ if(_me==null) _me=this;
+ _EMsampWt = 1.0/0.022;
+ _HDsampWt = 1.0/0.03;
+ _EMmip = 0.008;
+ _HDmip = 0.031;
+ _EMthresh = 0.25;
+ _HDthresh = 0.25;
+ _lyrNeighEM = 8; // was 8 ???
+ _zNeighEM = 4; // was 4 ???
+ _phiNeighEM = 4; // was 4 ???
+ _lyrNeighHD = 8; // was 8 ???
+ _zNeighHD = 4; // was 4 ???
+ _phiNeighHD = 4; // was 4 ???
+ _lyrContractionEM = new int[3];
+ _lyrContractionEM[0] = 0;
+ _lyrContractionEM[1] = -7;
+ _lyrContractionEM[2] = -7;
+ _zContractionEM = new int[3];
+ _zContractionEM[0] = -4;
+ _zContractionEM[1] = 2;
+ _zContractionEM[2] = 2;
+ _phiContractionEM = new int[3];
+ _phiContractionEM[0] = -6;
+ _phiContractionEM[1] = 2;
+ _phiContractionEM[2] = 2;
+ _lyrContractionHD = new int[3];
+ _zContractionHD = new int[3];
+ _phiContractionHD = new int[3];
+ _distanceType = "Euclidean";
+ _clusterSeparately = true;
+ _centroidWeightType = "Density";
+ _ModeValleyFactors = new double[2][2];
+ _ModeValleyFactors[0][0] = 0.8;
+ _ModeValleyFactors[0][1] = 0.15;
+ _ModeValleyFactors[1][0] = 0.025;
+ _ModeValleyFactors[1][1] = 0.025;
+ _convergenceParameter = 0.001;
+ _prune = true;
+ _pruningDist = 500.;
+ _maxMaskSize = 20;
+ _emSigma = new double[3];
+ _emSigma[0] = 1.5;
+ _emSigma[1] = 1.5;
+ _emSigma[2] = 1.5;
+ }
+
+ public double getEMweight(){
+ return this._EMsampWt;
+ }
+
+ public double getEMmip(){
+ return this._EMmip;
+ }
+
+ public double getEMthresh(){
+ return this._EMthresh;
+ }
+
+ public double getHDweight(){
+ return this._HDsampWt;
+ }
+
+ public double getHDmip(){
+ return this._HDmip;
+ }
+
+ public double getHDthresh(){
+ return this._HDthresh;
+ }
+
+ public int getLyrNeighEM(){
+ return this._lyrNeighEM;
+ }
+
+ public int getZNeighEM(){
+ return this._zNeighEM;
+ }
+
+ public int getPhiNeighEM(){
+ return this._phiNeighEM;
+ }
+
+ public int getLyrNeighHD(){
+ return this._lyrNeighHD;
+ }
+
+ public int getZNeighHD(){
+ return this._zNeighHD;
+ }
+
+ public int getPhiNeighHD(){
+ return this._phiNeighHD;
+ }
+
+ public String getDistanceType(){
+ return this._distanceType;
+ }
+
+ public String getCentroidWeightType(){
+ return this._centroidWeightType;
+ }
+
+ public boolean ClusterSeparately(){
+ return this._clusterSeparately;
+ }
+
+ public double[][] getInfluenceFactors(){
+ return this._ModeValleyFactors;
+ }
+
+ public double getConvergenceParameter(){
+ return this._convergenceParameter;
+ }
+
+ public boolean Prune(){
+ return this._prune;
+ }
+
+ public double getPruningDist(){
+ return this._pruningDist;
+ }
+
+ public int getMaxMaskSize(){
+ return this._maxMaskSize;
+ }
+
+ public double[] getEMGaussWidths(){
+ return this._emSigma;
+ }
+
+ public int[] getLyrContracEM(){
+ return this._lyrContractionEM;
+ }
+
+ public int[] getZContracEM(){
+ return this._zContractionEM;
+ }
+
+ public int[] getPhiContracEM(){
+ return this._phiContractionEM;
+ }
+
+ public int[] getLyrContracHD(){
+ return this._lyrContractionHD;
+ }
+
+ public int[] getZContracHD(){
+ return this._zContractionHD;
+ }
+
+ public int[] getPhiContracHD(){
+ return this._phiContractionHD;
+ }
+
+ private static RunControlParameters _me=null;
+
+ private double _EMsampWt;
+ private double _HDsampWt;
+ private double _EMmip;
+ private double _HDmip;
+ private double _EMthresh;
+ private double _HDthresh;
+ private double _convergenceParameter;
+ private double _pruningDist;
+ private int _lyrNeighEM,_lyrNeighHD;
+ private int _zNeighEM,_zNeighHD;
+ private int _phiNeighEM,_phiNeighHD;
+ private int _maxMaskSize;
+ private String _distanceType;
+ private String _centroidWeightType;
+ private boolean _clusterSeparately;
+ private boolean _prune;
+ private double[][] _ModeValleyFactors;
+ private double[] _emSigma;
+ private int[] _lyrContractionEM,_lyrContractionHD;
+ private int[] _zContractionEM,_zContractionHD;
+ private int[] _phiContractionEM,_phiContractionHD;
+
+}