13 added files
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N CalculateDistance.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CalculateDistance.java 8 Jan 2006 14:28:21 -0000 1.1
@@ -0,0 +1,58 @@
+package org.lcsim.recon.cluster.directedtree;
+//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;
+ }
+
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N Cluster.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Cluster.java 8 Jan 2006 14:28:21 -0000 1.1
@@ -0,0 +1,22 @@
+package org.lcsim.recon.cluster.directedtree;
+
+import java.util.Vector;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+public class Cluster extends BasicCluster {
+
+ public Cluster(double ene,double xpos,double ypos,double zpos,
+ Vector<CalorimeterHit> hitColl)
+ {
+ _ene = ene;
+ _xpos = xpos;
+ _ypos = ypos;
+ _zpos = zpos;
+ for( CalorimeterHit hit : hitColl ) {
+ addHit( hit );
+ }
+ }
+
+ double _ene,_xpos,_ypos,_zpos;
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N ClusterBuilder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClusterBuilder.java 8 Jan 2006 14:28:21 -0000 1.1
@@ -0,0 +1,35 @@
+package org.lcsim.recon.cluster.directedtree;
+import java.util.Vector;
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+public class ClusterBuilder {
+
+ HitWeightingClusterPropertyCalculator _hitWeightingCPC = new HitWeightingClusterPropertyCalculator();
+
+ public List<BasicCluster> makeClusters(Vector<Vector<CalorimeterHit>> inClusters)
+ {
+ int nclus = inClusters.size();
+ assert nclus > 0 : "ClusterBuilder: no clusters received.";
+
+ List<BasicCluster> outClusters = new ArrayList<BasicCluster>();
+
+ // loop over clusters
+ for( Vector<CalorimeterHit> hits : inClusters ){
+ int cluSize = hits.size();
+ assert cluSize > 0 : "ClusterBuilder: zero-hits cluster found.";
+
+ // create cluster
+ BasicCluster clus = new BasicCluster();
+ for( CalorimeterHit hit : hits ) clus.addHit( hit );
+
+ // change cluster's hit-weighting scheme
+ clus.setPropertyCalculator( _hitWeightingCPC );
+
+ outClusters.add(clus);
+ }
+ return outClusters;
+ }
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N DTreeAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DTreeAnalysis.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,449 @@
+package org.lcsim.recon.cluster.directedtree;
+
+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.digisim.CellSelector;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.SortClustersBySize;
+import org.lcsim.geometry.IDDecoder;
+
+public class DTreeAnalysis extends Driver {
+
+ private int _debug=1;
+
+ public DTreeAnalysis() {
+ _runPar = RunControlParameters.getInstance();
+ _emEcut = _runPar.getEMmip() * _runPar.getEMthresh() / _runPar.getEMweight();
+ _hdEcut = _runPar.getHDmip() * _runPar.getHDthresh() / _runPar.getHDweight();
+
+ // aida stuff
+ AIDA aida = AIDA.defaultInstance();
+ _tree = aida.tree();
+ _tf = aida.analysisFactory().createTupleFactory(_tree);
+ _tree.mkdir("DirectedTree");
+ _tree.cd("DirectedTree");
+ bookDirectedTreeNtuple();
+ _tree.cd("..");
+ }
+
+ public void process(EventHeader event) {
+
+ _evtnum = event.getEventNumber();
+ if(_debug>0) {
+ System.out.println("Starting DTreeAnalysis, event # "+_evtnum);
+ }
+ _tree.cd("DirectedTree");
+
+ String embName = "EcalBarrHits";
+ String emeName = "EcalEndcapHits";
+ String hdbName = "HcalBarrHits";
+ String hdeName = "HcalEndcapHits";
+ embhitmap = _expert.getCollHitMap(embName, _emEcut);
+ emehitmap = _expert.getCollHitMap(emeName, _emEcut);
+ hdbhitmap = _expert.getCollHitMap(hdbName, _hdEcut);
+ hdehitmap = _expert.getCollHitMap(hdeName, _hdEcut);
+
+ if(_debug>0) {
+ System.out.println("DTree: #hits: "
+ +" EMB="+embhitmap.size()
+ +", EMEC="+emehitmap.size()
+ +", HDB="+hdbhitmap.size()
+ +", HDEC="+hdehitmap.size());
+ }
+
+ if(_runPar.ClusterSeparately()){
+ _calType = "EM";
+ _nClusEM = 0; _nMCClusEM = 0;
+ processCollection(event,embName,embhitmap);
+ processCollection(event,emeName,emehitmap);
+
+ _calType = "HD";
+ _nClusHD = 0; _nMCClusHD = 0;
+ processCollection(event,hdbName,hdbhitmap);
+ processCollection(event,hdeName,hdehitmap);
+
+ finalizeTuple();
+ }
+ 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,int ily,int iz,int iphi}","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,int ily,int iz,int iphi}","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,float vx,float vy,float vz}","emmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi,float r,float phi,float theta,float time}","evtno"};
+ Class[] 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,float chrg,float e,float px,float py,float pz,float vx,float vy,float vz}","hdmcCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTagG,int cellTagR,int ily,int iz,int iphi,float r,float phi,float theta,float time}","evtno"};
+ Class[] columnClassesMChd = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
+ _tupleMChd = _tf.create("MCHD","Gen Clusters in HD",columnNamesMChd,columnClassesMChd);
+ }
+ else{
+ }
+ }
+
+ public void fillDirectedTreeNtuple(List<BasicCluster> recon,
+ List<BasicCluster> mctruth,
+ IDDecoder decoder)
+ {
+ 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);
+ }
+
+ // Fill tuples with reco info
+ int tag = 0;
+ if(_calType=="EM") tag = _nClusEM;
+ if(_calType=="HD") tag = _nClusHD;
+ Map<CalorimeterHit,Integer> recTagMap
+ = new HashMap<CalorimeterHit,Integer>();
+ for(BasicCluster clust : recon ) {
+ List<CalorimeterHit> cellVec = clust.getCalorimeterHits();
+// System.out.println("reco tag: "+tag+", #hits="+cellVec.size());
+ if(cellVec.size()==0)
+ System.out.println("*** bad cellVec in "+_calType);
+ for(CalorimeterHit ihit : cellVec) {
+ Integer tagObject = new Integer(tag);
+// System.out.println("putting <"+MyTools.printID(ihit.getCellID())+" -> "+tagObject+">");
+ recTagMap.put(ihit,tagObject);
+ double[] pos = ihit.getPosition();
+ double cdens = _loader.getDensity(ihit);
+ long jid = ihit.getCellID();
+ decoder.setID( ihit.getCellID() );
+ int ily = decoder.getValue(0); // layer
+ int iu = decoder.getValue(3); // either theta, phi or x
+ int iv = decoder.getValue(4); // either phi, z or y
+
+// int ily1 = MyTools.getLayer(jid);
+// int iz = MyTools.getThetaBin(jid);
+// int iphi = MyTools.getPhiBin(jid);
+// System.out.println("ilay="+ily1+" "+ily+", iu="+iphi+" "+iu+", iv="+iz+" "+iv);
+// assert ily==ily1 : "problem in layer index";
+// assert iphi==iu : "problem in phi/iU index";
+// assert iz==iv : "problem in iz/iV index";
+
+ ITuple cellFolder = null;
+ if(_calType=="EM") cellFolder = emCellFolder;
+ if(_calType=="HD") cellFolder = hdCellFolder;
+ if(cellFolder!=null) {
+ cellFolder.fill(0,(float) ihit.getRawEnergy());
+ cellFolder.fill(1,(float)pos[0]);
+ cellFolder.fill(2,(float)pos[1]);
+ cellFolder.fill(3,(float)pos[2]);
+ cellFolder.fill(4,(float)cdens);
+ cellFolder.fill(5,tag);
+ cellFolder.fill(6,ily);
+ cellFolder.fill(7,iu);
+ cellFolder.fill(8,iv);
+ cellFolder.addRow();
+ }
+ }
+ ITuple clusFolder = null;
+ if(_calType=="EM") clusFolder = emClusFolder;
+ if(_calType=="HD") clusFolder = hdClusFolder;
+ if(clusFolder!=null) {
+ Hep3Vector pos = new BasicHep3Vector( clust.getPosition() );
+ clusFolder.fill(0,(float) clust.getRawEnergy());
+ clusFolder.fill(1,(float) pos.x());
+ clusFolder.fill(2,(float) pos.y());
+ clusFolder.fill(3,(float) pos.z());
+ clusFolder.fill(4, clust.getSize());
+ clusFolder.addRow();
+ }
+ tag++;
+ }
+ if(_calType=="EM") _nClusEM = tag;
+ if(_calType=="HD") _nClusHD = tag;
+// System.out.println("calType="+_calType+", tag="+tag+", nrecoEM="+_nClusEM+", nrecoHD="+_nClusHD);
+
+ // Fill tuples with MC info
+ tag = 0;
+ if(_calType=="EM") tag = _nMCClusEM;
+ if(_calType=="HD") tag = _nMCClusHD;
+ for( BasicCluster clust : mctruth ) {
+ for( CalorimeterHit ihit : clust.getCalorimeterHits() ) {
+ Integer tagObject = recTagMap.get(ihit);
+
+ int rtag = -1;
+ if(tagObject!=null) rtag = tagObject.intValue();
+// int rtag = tagObject.intValue();
+// System.out.println("getting <"+MyTools.printID(ihit.getCellID())+" -> "+rtag+">");
+ double[] pos = ihit.getPosition();
+ double cdens = _loader.getDensity(ihit);
+ if(cdens==0) continue; // skip hits below energy cut
+ long jid = ihit.getCellID();
+ decoder.setID( ihit.getCellID() );
+ int ily = decoder.getValue(0); // layer
+ int iu = decoder.getValue(3); // either theta, phi or x
+ int iv = decoder.getValue(4); // either phi, z or y
+
+// int ily1 = MyTools.getLayer(jid);
+// int iz = MyTools.getThetaBin(jid);
+// int iphi = MyTools.getPhiBin(jid);
+// System.out.println("ilay="+ily1+" "+ily+", iu="+iz+" "+iu+", iv="+iphi+" "+iv);
+// assert ily==ily1 : "problem in layer index";
+// assert iphi==iv : "problem in phi/iV index";
+// assert iz==iu : "problem in iz/iU index";
+
+ double cellR = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
+ double theta = Math.atan2(cellR,pos[2]);
+ double phi = Math.atan2(pos[1],pos[0]);
+ if(phi<0.) phi += 2.0*Math.PI;
+ SimCalorimeterHit iihit = (SimCalorimeterHit)ihit;
+ int npart = iihit.getMCParticleCount();
+ double timeStamp = 99999.;
+ for(int j=0;j<npart;j++) {
+ double tmpTime = (double) iihit.getContributedTime(j);
+ if(tmpTime<timeStamp){
+ timeStamp = tmpTime;
+ }
+ }
+
+ ITuple mcCellFolder = null;
+ if(_calType=="EM") mcCellFolder = emmcCellFolder;
+ if(_calType=="HD") mcCellFolder = hdmcCellFolder;
+ if(mcCellFolder!=null) {
+ mcCellFolder.fill(0,(float) ihit.getRawEnergy());
+ mcCellFolder.fill(1,(float)pos[0]);
+ mcCellFolder.fill(2,(float)pos[1]);
+ mcCellFolder.fill(3,(float)pos[2]);
+ mcCellFolder.fill(4,(float) cdens);
+// System.out.println("Gtag="+tag+", Rtag="+rtag);
+ mcCellFolder.fill(5,tag);
+ mcCellFolder.fill(6,rtag);
+ mcCellFolder.fill(7,ily);
+ mcCellFolder.fill(8,iu);
+ mcCellFolder.fill(9,iv);
+ mcCellFolder.fill(10,(float) cellR);
+ mcCellFolder.fill(11,(float) phi);
+ mcCellFolder.fill(12,(float) theta);
+ mcCellFolder.fill(13,(float) timeStamp);
+ mcCellFolder.addRow();
+// System.out.println("EM Clus cells tag="+tag+", rtag="+rtag
+// +", id="+MyTools.printID(jid)
+// // +", pos="+pos[0]+" "+pos[1]+" "+pos[2]
+// +", dens="+cdens);
+ }
+ }
+
+ MCParticle mcpart = getMCParticleInCluster(clust);
+ if(mcpart!=null) {
+ float charge = (float)mcpart.getCharge();
+ int particleID = mcpart.getType().getPDGID();
+ float mcE = (float)mcpart.getEnergy();
+ Hep3Vector mcMom = mcpart.getMomentum();
+ double[] pos = clust.getPosition();
+ Hep3Vector vert = mcpart.getOrigin();
+
+ ITuple mcClusFolder = null;
+ if(_calType=="EM") mcClusFolder = emmcClusFolder;
+ if(_calType=="HD") mcClusFolder = hdmcClusFolder;
+ if(mcClusFolder!=null) {
+ mcClusFolder.fill(0,(float) clust.getRawEnergy());
+ mcClusFolder.fill(1, (float)pos[0] );
+ mcClusFolder.fill(2, (float)pos[1] );
+ mcClusFolder.fill(3, (float)pos[2] );
+ mcClusFolder.fill(4, clust.getSize() );
+ mcClusFolder.fill(5, particleID);
+ mcClusFolder.fill(6, charge);
+ mcClusFolder.fill(7, mcE);
+ mcClusFolder.fill(8, (float)mcMom.x());
+ mcClusFolder.fill(9,(float)mcMom.y());
+ mcClusFolder.fill(10,(float)mcMom.z());
+ mcClusFolder.fill(11,(float) vert.x());
+ mcClusFolder.fill(12,(float) vert.y());
+ mcClusFolder.fill(13,(float) vert.z());
+ mcClusFolder.addRow();
+ }
+ tag++;
+ }
+ }
+ if(_calType=="EM") _nMCClusEM = tag;
+ if(_calType=="HD") _nMCClusHD = tag;
+// System.out.println("calType="+_calType+", tag="+tag+", nMCEM="+_nMCClusEM+", nMCHD="+_nMCClusHD);
+
+ recTagMap.clear();
+ }
+
+ public void finalizeTuple() {
+ // EM
+ _tupleEM.fill(0,_nClusEM);
+ _tupleEM.fill(3,_evtnum);
+ _tupleEM.addRow();
+ _tupleMCem.fill(0,_nMCClusEM);
+ _tupleMCem.fill(3,_evtnum);
+ _tupleMCem.addRow();
+ // HD
+ _tupleHD.fill(0,_nClusHD);
+ _tupleHD.fill(3,_evtnum);
+ _tupleHD.addRow();
+ _tupleMChd.fill(0,_nMCClusHD);
+ _tupleMChd.fill(3,_evtnum);
+ _tupleMChd.addRow();
+ // event
+ _tupleEvt.fill(0,_evtnum);
+ _tupleEvt.addRow();
+ }
+
+ private void processCollection(EventHeader event, String colName,
+ Map<Long,CalorimeterHit> hitmap)
+ {
+ if(hitmap.size()==0) return;
+ IDDecoder decoder = _expert.getIDDecoder(colName);
+
+ // get NearestNeighbor clusters
+ List<BasicCluster> nnClusColl = new ArrayList<BasicCluster>();
+ try {
+ nnClusColl = event.get(BasicCluster.class, colName+"NNClusters");
+ }
+ catch(IllegalArgumentException x) {
+ // ignore
+ }
+
+ // get reco clusters
+ List<BasicCluster> recoClusColl = new ArrayList<BasicCluster>();
+ 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() );
+ if(_debug>0) {
+ System.out.println("*** DTreeAnalysis: colName="+colName
+ +", # of MCclusters: "+mcClusColl.size()
+ +", # of NNclusters: "+nnClusColl.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();
+
+ // Analyze NN clusters
+ int nnhits=0;
+ for( BasicCluster iclus : nnClusColl ) {
+ double pos[] = iclus.getPosition();
+ List<CalorimeterHit> iHits = iclus.getCalorimeterHits();
+ nnhits += iHits.size(); // includes hits below threshold
+ }
+ int nNNclus = nnClusColl==null ? 0 : nnClusColl.size();
+
+ // Analyze DirectedTree clusters
+ int dthits=0;
+ for( BasicCluster iclus : recoClusColl ) {
+ double pos[] = iclus.getPosition();
+ List<CalorimeterHit> iHits = iclus.getCalorimeterHits();
+ dthits += iHits.size(); // includes hits below threshold
+ }
+ int nDTclus = recoClusColl==null ? 0 : recoClusColl.size();
+
+ if(_debug>0) {
+ System.out.println(" Comparing #clusters: "
+ +" MCclus="+nclusCheat
+ +" NNclus="+nNNclus
+ +" DTclus="+nDTclus
+ );
+ System.out.println(" Comparing total #hits: "
+ +" in event: "+hitmap.size()
+ +", in MCclus="+mchits
+ +", in NNclus="+nnhits
+ +", in DTclus="+dthits );
+ }
+
+ 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, decoder );
+ }
+
+ 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 CalHitMapMgr _expert = CalHitMapMgr.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();
+ private int layerIndex, uIndex, vIndex;
+ private double _emEcut, _hdEcut;
+
+ private int _nClusEM,_nMCClusEM,_nClusHD,_nMCClusHD;
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N DirectedTreeClusterer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DirectedTreeClusterer.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,398 @@
+package org.lcsim.recon.cluster.directedtree;
+
+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.CalHitMapMgr;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.SortClustersBySize;
+import org.lcsim.event.Cluster;
+
+public class DirectedTreeClusterer extends Driver {
+
+ private int _debug = 0;
+
+ public DirectedTreeClusterer() {
+ _rcp = RunControlParameters.getInstance();
+ _emEcut = _rcp.getEMmip() * _rcp.getEMthresh() / _rcp.getEMweight();
+ _hdEcut = _rcp.getHDmip() * _rcp.getHDthresh() / _rcp.getHDweight();
+ }
+
+ public void process(EventHeader event) {
+
+ _evtnum = event.getEventNumber();
+ if(_debug>0) {
+ System.out.println("Start DirectedTreeClusterer, event #"+_evtnum);
+ }
+
+ _roMap = event.getDetector().getReadouts();
+
+ String embName = "EcalBarrHits";
+ String emeName = "EcalEndcapHits";
+ String hdbName = "HcalBarrHits";
+ String hdeName = "HcalEndcapHits";
+ embhitmap = _expert.getCollHitMap(embName, _emEcut);
+ emehitmap = _expert.getCollHitMap(emeName, _emEcut);
+ hdbhitmap = _expert.getCollHitMap(hdbName, _hdEcut);
+ hdehitmap = _expert.getCollHitMap(hdeName, _hdEcut);
+ _loader.setDensities(embName, embhitmap);
+ _loader.setDensities(emeName, emehitmap);
+ _loader.setDensities(hdbName, hdbhitmap);
+ _loader.setDensities(hdeName, hdehitmap);
+
+ if(_debug>0) {
+ System.out.println("DTree: #hits: EMB="+embhitmap.size()
+ +", EMEC="+emehitmap.size()
+ +", HB="+hdbhitmap.size()
+ +", HEC="+hdehitmap.size());
+ }
+
+ if(_rcp.ClusterSeparately()) {
+ _calType = "EM";
+ processCollection(event, embName, embhitmap);
+ processCollection(event, emeName, emehitmap);
+
+ _calType = "HD";
+ processCollection(event, hdbName, hdbhitmap);
+ processCollection(event, hdeName, hdehitmap);
+ }
+ 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,_rcp);
+ }
+
+ }
+
+ public Vector<Vector<CalorimeterHit>>
+ makeTree(Map<Long,CalorimeterHit> cahitmap, SegmentationBase segm) {
+
+ if(_debug>0) {
+ System.out.println("makeTree: 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>();
+ // changed from !=0
+ if(idens>=0.0) {
+ double[] ipos = ihit.getPosition();
+ if(_calType=="EM"){
+ int nLyrOrig = _rcp.getLyrNeighEM();
+ int nZOrig = _rcp.getZNeighEM();
+ int nPhiOrig = _rcp.getPhiNeighEM();
+ int dRegion = emDensityRegion(idens,nLyrOrig,nZOrig,nPhiOrig);
+ int[] lyrCon = _rcp.getLyrContracEM();
+ int[] zCon = _rcp.getZContracEM();
+ int[] phiCon = _rcp.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 = _rcp.getLyrNeighHD();
+ int nZOrig = _rcp.getZNeighHD();
+ int nPhiOrig = _rcp.getPhiNeighHD();
+ int dRegion = hdDensityRegion(idens,nLyrOrig,nZOrig,nPhiOrig);
+ int[] lyrCon = _rcp.getLyrContracHD();
+ int[] zCon = _rcp.getZContracHD();
+ int[] phiCon = _rcp.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);
+
+ 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 = _rcp.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(ihit);
+// 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>();
+ 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) {
+ CalorimeterHit jhit = cahitmap.get(jkey);
+ double[] jpos = jhit.getPosition();
+ CalculateDistance calcD = new CalculateDistance();
+ double distance = calcD.CalculateDistance(_distType,ipos,jpos);
+ if(distance<dmin){
+ dmin = distance;
+ parentKey = jkey;
+ }
+ }
+ _ParentMap.put(cellid,parentKey);
+// System.out.println("Cell is matched to identical density neighbor");
+ }
+ }
+ }
+ else{
+ RootVec.add(ihit);
+// System.out.println("Zero density root");
+ }
+ }
+
+ 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);
+ }
+ }
+ if(_debug>0) {
+ System.out.println("# starting points = "+startingPoints.size());
+ }
+
+ Vector<Vector<CalorimeterHit>> branches = new Vector<Vector<CalorimeterHit>>();
+ for(int i=0;i<startingPoints.size();i++) {
+ 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();
+ if(_debug>0) System.out.println("no. of roots = "+vsiz);
+ Vector<Vector<CalorimeterHit>> caTrees = new Vector<Vector<CalorimeterHit>>();
+ for(int i=0;i<vsiz;i++){
+ caTrees.add( new Vector<CalorimeterHit>() );
+ }
+
+ 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(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(_debug>0) {
+ if(caTrees.get(i).size()>1) {
+ System.out.println("ROOT id="+MyTools.printID(tmpID)
+ +", #hits="+caTrees.get(i).size());
+ }
+ else ++nTrivial; // for single-hit clusters
+ }
+ }
+ if(_debug>0) {
+ 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<6.) region=1;
+ if(dens>=6. && dens<26.) region=2;
+ if(dens>=26.) region=3;
+ // }
+ return region;
+ }
+
+ public int hdDensityRegion(double dens,int nl,int nz,int np){
+ int region = 0;
+ if(dens<=6.) region=1;
+ if(dens>=6. && dens<26.) region=2;
+ if(dens>=26.) region=3;
+ return region;
+ }
+
+ private void processCollection(EventHeader event, String colName,
+ Map<Long,CalorimeterHit> hitmap) {
+
+ if(_debug>0) {
+ 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");
+ event.put( newName, recoClusColl, Cluster.class, (1<<31) );
+ }
+
+// 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 _rcp;
+ private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+ private CalHitMapMgr _expert = CalHitMapMgr.getInstance();
+ private Map<Long,Long> _ParentMap;
+ private Map<Long,CalorimeterHit> embhitmap,emehitmap,hdbhitmap,hdehitmap;
+ private String _distType,_calType;
+ private int _nLyr;
+ private int _nZ;
+ private int _nPhi;
+ private int _evtnum;
+ private double _emEcut, _hdEcut;
+
+ private ClusterBuilder clusBuilder = new ClusterBuilder();
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N HitWeightingClusterPropertyCalculator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HitWeightingClusterPropertyCalculator.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,84 @@
+package org.lcsim.recon.cluster.directedtree;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.recon.cluster.util.DefaultClusterPropertyCalculator;
+
+/**
+ * Density weighting implementation of a ClusterPropertyCalculator object.
+ *
+ * <br>This alternative hit-weighting scheme can be used for new cluster
+ * properties in two significantly different ways:
+ *
+ * <ol>
+ * <li> By permanently changing the weighting scheme of the Cluster objects:
+ * <pre>
+ * aCluster.setPropertyCalculator( myCalculator );
+ * Hep3Vector newPos = new BasicHep3Vector( aCluster.getPosition() );
+ * </pre></li>
+ *
+ * <li> Using another weighting scheme without altering the scheme used
+ * by the Cluster objects:
+ * <pre>
+ * myCalculator.calculateProperties( cluster.getCalorimeterHits() );
+ * Hep3Vector newPos = new BasicHep3Vector( myCalculator.getPosition() );
+ * </pre></li>
+ * </ol>
+ *
+ * @author Guilherme Lima
+ */
+public class HitWeightingClusterPropertyCalculator extends DefaultClusterPropertyCalculator
+{
+ RunControlParameters _runPar = RunControlParameters.getInstance();
+ LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+
+ public void calculateProperties(List<CalorimeterHit> hits)
+ {
+ double cluE = 0.0;
+ double sumWt = 0.0;
+ double cluX = 0.0; double cluY = 0.0; double cluZ = 0.0;
+ int cluSize = hits.size();
+ if(cluSize==0) System.out.println("Damn 2");
+ assert cluSize > 0 : "ClusterBuilder: zero-hits cluster found.";
+
+ for( CalorimeterHit hit : hits ){
+ double ene = hit.getRawEnergy();
+ double[] pos = hit.getPosition();
+ cluE += ene;
+ if(_runPar.getCentroidWeightType()=="Energy") {
+ cluX += ene*pos[0];
+ cluY += ene*pos[1];
+ cluZ += ene*pos[2];
+ sumWt += ene;
+ }
+ if(_runPar.getCentroidWeightType()=="Density") {
+// if(cluSize>1) {
+ long id = hit.getCellID();
+ double dens = _loader.getDensity(hit.getCellID());
+ 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 /= sumWt;
+ cluY /= sumWt;
+ cluZ /= sumWt;
+ }
+
+ position[0] = cluX;
+ position[1] = cluY;
+ position[2] = cluZ;
+ }
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N Kinem.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Kinem.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,18 @@
+package org.lcsim.recon.cluster.directedtree;
+
+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 += 2.0*Math.PI;
+ return phi;
+ }
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N LoadMyCalorimeterHit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LoadMyCalorimeterHit.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,110 @@
+package org.lcsim.recon.cluster.directedtree;
+
+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.CalorimeterHit;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.segmentation.SegmentationBase;
+import org.lcsim.geometry.compact.Readout;
+import org.lcsim.geometry.Subdetector;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
+import org.lcsim.geometry.IDDecoder;
+
+public class LoadMyCalorimeterHit {
+
+ public static LoadMyCalorimeterHit getInstance() {
+ if(_me==null) {
+ new LoadMyCalorimeterHit();
+ System.out.println("LoadMyCalorimeterHit object at "+_me);
+ }
+ assert _me != null : "Problem constructing LoadMyCalorimeterHit.";
+ return _me;
+ }
+
+ private LoadMyCalorimeterHit() {
+ if(_me==null) _me=this;
+ _densityMap = new HashMap<Long,Double>();
+ _runPar = RunControlParameters.getInstance();
+ }
+
+ public void setEvent(EventHeader event) {
+ if( event != _event ) {
+ this.reset();
+ }
+ }
+
+ public void reset() {
+ // clear hit maps for each subdetector collection
+ _densityMap.clear();
+ }
+
+ // assign densities to every hit in *new* hitmap,
+ // thus density=0 for hits failing energy cut
+ public void setDensities(String colName, Map<Long,CalorimeterHit> hitmap) {
+
+ if( colName.contains("Ecal") ) {
+ _thresh = _runPar.getEMthresh();
+ _nLyr = _runPar.getLyrNeighEM();
+ _nZ = _runPar.getZNeighEM();
+ _nPhi = _runPar.getPhiNeighEM();
+ }
+ else if( colName.contains("Hcal") ) {
+ _thresh = _runPar.getHDthresh();
+ _nLyr = _runPar.getLyrNeighHD();
+ _nZ = _runPar.getZNeighHD();
+ _nPhi = _runPar.getPhiNeighHD();
+ }
+
+ for( CalorimeterHit khit : hitmap.values() ) {
+ long cellid = khit.getCellID();
+ IDDecoder decoder = khit.getSubdetector().getIDDecoder();
+ SegmentationBase segm = (SegmentationBase)decoder;
+
+ double nNeigh = 1.0; // count hit itself + neighbors (never zero)
+ segm.setID(cellid);
+ long neigh[] = segm.getNeighbourIDs(_nLyr,_nZ,_nPhi);
+ for(int j=0; j<neigh.length; ++j) {
+ long jid = neigh[j];
+ CalorimeterHit lhit = hitmap.get(jid);
+ if(lhit!=null) nNeigh++;
+ }
+
+ _densityMap.put(cellid,nNeigh);
+ if(nNeigh<50) AIDA.defaultInstance().cloud1D("density").fill(nNeigh);
+ }
+ }
+
+ // Get hit density from a hit reference
+ double getDensity(final CalorimeterHit hit) {
+ if(hit==null) return 0.0;
+ return this.getDensity( hit.getCellID() );
+ }
+
+ // Get hit density from cellID
+ double getDensity(long cellid) {
+ try {
+ return _densityMap.get( cellid );
+ }
+ catch(NullPointerException e) {
+ // Hits below threshold end up here, should return density=0
+ return 0.0;
+ }
+ }
+
+ //*** FIELDS ***
+
+ private static LoadMyCalorimeterHit _me = null;
+ /** Current event */
+ private EventHeader _event;
+ private static RunControlParameters _runPar;
+ private double _thresh;
+ private int _nLyr;
+ private int _nZ;
+ private int _nPhi;
+ private Map<Long,Double> _densityMap;
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N MyTools.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MyTools.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,117 @@
+/*
+ * MyTools.java
+ *
+ * Created on Nov. 2003
+ * G.Lima, J. McCormick (NICADD/NIU)
+ */
+
+package org.lcsim.recon.cluster.directedtree;
+import org.lcsim.event.EventHeader;
+
+/**
+ *
+ * @author G.Lima, J.McCormick
+ */
+public class MyTools
+{
+ // 1=sidaug05; 2=sidaug05_tcmt; 3=sdnphoct04
+ static int detector = 2;
+
+// // 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 getSystem( long cellid ) {
+ if(detector==1) return (int)(cellid>>7)&0x3f;
+ else if(detector==2) return (int)(cellid>>7)&0x3f;
+ else if(detector==3) return (int)(cellid>>16)&0x3f;
+ else return 0;
+ }
+
+ public static int getSystemBarrel( long cellid ) {
+ if(detector==1) return (int)((cellid>>7) & 0x1ff);
+ else if(detector==2) return (int)((cellid>>7) & 0x1ff);
+ else if(detector==3) return (int)((cellid>>16) & 0x1ff);
+ else return 0;
+ }
+
+ public static int getLayer( long cellid ) {
+ if(detector==1) return (int)((cellid>>0) & 0x7f);
+ else if(detector==2) return (int)((cellid>>0) & 0x7f);
+ else if(detector==3) return (int)((cellid>>0) & 0xffff);
+ else return 0;
+ }
+
+ public static int getThetaBin( long cellid ) {
+ int thetabin = 0;
+ if(detector==1) {
+ int max = 1<<11; // for 11 bits
+ thetabin = (int)((cellid>>32) & 0x7ff);
+ }
+ if(detector==2) {
+ int max = 1<<16; // for 16 bits
+ thetabin = (int)((cellid>>48) & 0xffff);
+ // sign bit
+ if( thetabin > max/2-1 ) thetabin = thetabin - max;
+ }
+ if(detector==3) {
+ int max = 1<<16; // for 16 bits
+ thetabin = (int)((cellid>>32) & 0xffff);
+ // sign bit
+ if( thetabin > max/2-1 ) thetabin = thetabin - max;
+ }
+ return thetabin;
+ }
+
+ public static int getPhiBin( long cellid ) {
+ int phibin = 0;
+ if(detector==1) {
+ int max = 1<<11; // for 11 bits
+ phibin = (int)((cellid>>43) & 0x7ff);
+ if( phibin > max/2-1 ) phibin = phibin - max;
+ }
+ if(detector==2) phibin = (int)((cellid>>32) & 0xffff);
+ if(detector==3) phibin = (int)((cellid>>48) & 0xffff);
+ return phibin;
+ }
+
+ 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( getSystemBarrel(id) )
+ +" "+Integer.toString( getLayer(id) )
+ +" "+Integer.toString( getThetaBin(id))
+ +" "+Integer.toString( getPhiBin(id) )
+ +">");
+ }
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N RunControlParameters.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RunControlParameters.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,191 @@
+package org.lcsim.recon.cluster.directedtree;
+
+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; // in GeV (corrected energy)
+ _HDmip = 0.031; // in GeV (corrected energy)
+ _EMthresh = 0.25; // was 0.25
+ _HDthresh = 0.25; // was 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] = 0;
+ _zContractionEM[2] = 2;
+ _phiContractionEM = new int[3];
+ _phiContractionEM[0] = -6;
+ _phiContractionEM[1] = 0;
+ _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;
+
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N TrackHitMatcher.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackHitMatcher.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,150 @@
+package org.lcsim.recon.cluster.directedtree;
+
+import java.util.Map;
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+
+// import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Track;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.layer.Layering;
+import org.lcsim.geometry.segmentation.SegmentationBase;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
+
+/**
+ * A helper class for track-hit matching in a single calorimeter component
+ *
+ * @author Guilherme Lima
+ * @version $Id: TrackHitMatcher.java,v 1.1 2006/01/08 14:28:22 lima Exp $
+ */
+public class TrackHitMatcher {
+
+ int _debug = 0; // debug level, 0 for no printout
+
+ // constructors
+ public TrackHitMatcher(String colName, double[] layers) {
+ this(colName, layers, 2, 2, 0);
+ }
+ public TrackHitMatcher(String colName, double[] layers, int dU, int dV) {
+ this(colName, layers, dU, dV, 0);
+ }
+
+ public TrackHitMatcher(String colName, double[] layers, int dU, int dV, int debug)
+ {
+ _colName = colName;
+ _layers = layers;
+ _debug = debug;
+ _subdet = _expert.getSubdetector( colName );
+ _segm = (SegmentationBase)_expert.getIDDecoder( colName );
+ _dU = dU;
+ _dV = dV;
+ }
+
+ public void findMatches(EventHeader event,
+ Map<Track,Vector<Hep3Vector>> intersects,
+ Map<Track,List<CalorimeterHit>> trkHitsMap)
+ {
+ if(!_init) initialize(event);
+ this.reset();
+
+ // Load hit collections
+ _hitmap = _expert.getCollHitMap(_colName);
+
+ // loop over tracks
+ for( Track trk : event.get(Track.class,"CombinedTracks") ) {
+ _trk = trk;
+ _trkHitsMap = trkHitsMap;
+ _trkHits = trkHitsMap.get( trk );
+ int nhitsBefore = (_trkHits!=null ? _trkHits.size() : 0);
+
+ Hep3Vector pvec = new BasicHep3Vector( trk.getMomentum() );
+ Vector<Hep3Vector> trkIntercepts = intersects.get(trk);
+ if( trkIntercepts == null) continue;
+
+ for( Hep3Vector pos : trkIntercepts ) {
+ // find cell containing track-cylinder intersection
+ long cellid = _segm.findCellContainingXYZ( pos );
+ if(cellid==0) continue;
+
+ // save matched hits
+ this.addHitToTrack( cellid );
+
+ _segm.setID( cellid );
+// double rhoCell = _segm.getDistanceToSensitive(i);
+// System.out.println("layer="+i+", pos="+pos+", cell pos=("
+// +_segm.getX()+"; "+_segm.getY()+"; "
+// +_segm.getZ()+")"
+// +", rhoSwim="+Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y())
+// +", rhoCell="+rhoCell+", #hits="+_hitmap.size());
+
+ // save matched hits in neighborhood
+ long[] neighs = _segm.getNeighbourIDs( 0, _dU, _dV );
+ for(int i=0; i<neighs.length; ++i) {
+ this.addHitToTrack( neighs[i] );
+ }
+ }
+
+ if(_debug>0) {
+ int nhitsAfter = (_trkHits!=null ? _trkHits.size() : 0);
+ System.out.println(_colName+" matcher:"
+ +" pvec=("+pvec.x()+"; "+pvec.y()+"; "+pvec.z()+")"
+ +", #intersects="+intersects.get(trk).size()
+ +", #matches="+(nhitsAfter-nhitsBefore) );
+ }
+ }
+ }
+
+ private void reset() {
+ }
+
+ private void initialize(EventHeader event) {
+ // layer info
+ CylindricalCalorimeter embSubdet = (CylindricalCalorimeter)_expert.getSubdetector(_colName);
+ Layering layers = embSubdet.getLayering();
+ int nlayers = layers.getLayerCount();
+ _layers = new double[nlayers];
+ for(int i=0; i<nlayers; ++i) {
+ _layers[i] = layers.getDistanceToLayerSensorMid(i);
+ }
+
+ _init = true;
+ }
+
+ // Checks whether hit is valid
+ private void addHitToTrack( long cellid ) {
+ CalorimeterHit ihit = _hitmap.get(cellid);
+ if(ihit!=null) {
+ if(_trkHits==null) {
+ _trkHits = new ArrayList<CalorimeterHit>();
+ _trkHitsMap.put( _trk, _trkHits );
+ }
+ if(_debug>1) System.out.println(" add "+MyTools.printID(cellid));
+ _trkHits.add( ihit );
+ }
+ }
+
+ //***** FIELDS ****
+
+ // track currently being processed
+ private boolean _init = false;
+ private String _colName;
+ private double[] _layers;
+ private Subdetector _subdet = null;
+ private SegmentationBase _segm = null;
+ private Map<Long,CalorimeterHit> _hitmap = null;
+// private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+ private CalHitMapMgr _expert = CalHitMapMgr.getInstance();
+ private int _dU, _dV;
+// private AIDA _aida = AIDA.defaultInstance();
+ // these are for loop optimization
+ private Track _trk = null;
+ private List<CalorimeterHit> _trkHits = null;
+ private Map<Track,List<CalorimeterHit>> _trkHitsMap = null;
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N TrackHitRelation.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackHitRelation.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,18 @@
+package org.lcsim.recon.cluster.directedtree;
+
+import java.util.List;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Track;
+
+public class TrackHitRelation {
+ TrackHitRelation(Track track, List<CalorimeterHit> hits) {
+ _track = track;
+ _hits = hits;
+ }
+
+ public Track getTrack() { return _track; }
+ public List<CalorimeterHit> getHits() { return _hits; }
+
+ private Track _track;
+ private List<CalorimeterHit> _hits;
+}
lcsim/src/org/lcsim/recon/cluster/directedtree
diff -N TrackMatchingDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackMatchingDriver.java 8 Jan 2006 14:28:22 -0000 1.1
@@ -0,0 +1,237 @@
+package org.lcsim.recon.cluster.directedtree;
+
+import java.util.Map;
+import java.util.List;
+import java.util.Vector;
+import java.util.HashMap;
+import java.util.ArrayList;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Track;
+import org.lcsim.geometry.util.CalorimeterIDDecoder;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.segmentation.BarrelCylinderSegmentationBase;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.layer.Layering;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.util.swim.HelixSwim;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
+
+/**
+ * A driver for track matching using the swimmer.
+ *
+ * @author Guilherme Lima
+ * @version $Id: TrackMatchingDriver.java,v 1.1 2006/01/08 14:28:22 lima Exp $
+ */
+public class TrackMatchingDriver extends Driver {
+
+ int _debug = 2; // debug level, 0 for no printout
+
+ public TrackMatchingDriver() {
+ }
+
+ public void process(EventHeader event) {
+
+ if(_debug>0) {
+ System.out.println("******** Track to calorimeter hit matching ***");
+ }
+
+ if(!_init) initialize(event);
+ this.reset();
+
+ Map<Track,Vector<Hep3Vector>> emInterceptsMap = new HashMap<Track,Vector<Hep3Vector>>();
+ Map<Track,Vector<Hep3Vector>> hadInterceptsMap = new HashMap<Track,Vector<Hep3Vector>>();
+
+ // loop over tracks
+ List<Track> recoTracks = event.get(Track.class, "CombinedTracks");
+ System.out.println("TrackingCheater: # parts="+recoTracks.size());
+ for( Track trk : recoTracks ) {
+ Hep3Vector vtx = new BasicHep3Vector( trk.getReferencePoint() );
+ Hep3Vector pvec = new BasicHep3Vector( trk.getMomentum() );
+ int icharge = trk.getCharge();
+ if(_debug>1) {
+ System.out.println("*** New track"
+ +", pvec=("+pvec.x()+"; "+pvec.y()+"; "+pvec.z()+")"
+ +", vtx=("+vtx.x()+"; "+vtx.y()+"; "+vtx.z()+")"
+ +", q="+icharge);
+ }
+
+ // setup the swimmer for this track
+ _swimmer.setTrack( pvec, vtx, icharge );
+ double sParm1 = _swimmer.getDistanceToCylinder( 1270, 3000 );
+ Hep3Vector pos1 = _swimmer.getPointAtDistance( sParm1 );
+ double rho1 = Math.sqrt(pos1.x()*pos1.x()+pos1.y()*pos1.y());
+ System.out.println("Swimmer@EMentrance: ("+pos1.x()+"; "+pos1.y()+"; "+pos1.z()+", rho="+rho1);
+
+ // Swim the particle to each sensitive layer, saving intersections
+ boolean looping = false;
+ Vector<Hep3Vector> emIntercepts = new Vector<Hep3Vector>();
+ for(int i=0; i<_layersEMB.length; ++i) {
+ if(looping) break;
+
+ double rcyl = _layersEMB[i];
+ double zcyl = _layersEME[i];
+// System.out.println("*** Swimming to layer "+i
+// +", rcyl="+rcyl+", zcyl="+zcyl);
+ double sParm = _swimmer.getDistanceToCylinder( rcyl, zcyl );
+ Hep3Vector pos = _swimmer.getPointAtDistance( sParm );
+
+ // correction for plug-like endcaps
+ double rho = Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y());
+ if( rcyl-rho>1.e-3 && rho>_rhoMinEM ) {
+// System.out.println("Plug-type correction: pos="+pos+", rho="+rho);
+ sParm = _swimmer.getDistanceToRadius( rcyl );
+ pos = _swimmer.getPointAtDistance( sParm );
+ rho = Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y());
+ if(rho!=rcyl) {
+ looping = true;
+ if(_debug>1) System.out.println("Track seems to be looping. Break.");
+ break;
+ }
+ }
+ if(_debug>1) System.out.println("Swimmer: layer="+i+", pos="+pos);
+ emIntercepts.add(i,pos);
+ }
+
+ Vector<Hep3Vector> hadIntercepts = new Vector<Hep3Vector>();
+ for(int i=0; i<_layersHDB.length; ++i) {
+ if(looping) break;
+
+ double rcyl = _layersHDB[i];
+ double zcyl = _layersHDE[i];
+// System.out.println("*** Swimming to layer "+i
+// +", rcyl="+rcyl+", zcyl="+zcyl);
+ double sParm = _swimmer.getDistanceToCylinder( rcyl, zcyl );
+ Hep3Vector pos = _swimmer.getPointAtDistance( sParm );
+
+ // correction for plug-like endcaps
+ double rho = Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y());
+ if( rcyl-rho>1.e-3 && rho>_rhoMinHAD ) {
+// System.out.println("Plug-type correction: pos="+pos+", rho="+rho);
+ sParm = _swimmer.getDistanceToRadius( rcyl );
+ pos = _swimmer.getPointAtDistance( sParm );
+ rho = Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y());
+ if(rho!=rcyl) {
+ looping = true;
+ if(_debug>1) System.out.println("Track seems to be looping. Break.");
+ break;
+ }
+ }
+ if(_debug>1) System.out.println("Swimmer: layer="+i+", pos="+pos);
+ hadIntercepts.add(i,pos);
+ }
+ if(emIntercepts.size()>0) emInterceptsMap.put( trk, emIntercepts );
+ if(hadIntercepts.size()>0) hadInterceptsMap.put( trk, hadIntercepts );
+ }
+
+ // call track-hit matchers
+ Map<Track,List<CalorimeterHit>> trkHitsMap
+ = new HashMap<Track,List<CalorimeterHit>>();
+
+ _embMatcher.findMatches( event, emInterceptsMap, trkHitsMap );
+ _emeMatcher.findMatches( event, emInterceptsMap, trkHitsMap );
+ _hdbMatcher.findMatches( event, hadInterceptsMap, trkHitsMap );
+ _hdeMatcher.findMatches( event, hadInterceptsMap, trkHitsMap );
+
+ // summary
+ int i=0;
+ List<TrackHitRelation> trackHitRels = new ArrayList<TrackHitRelation>();
+ for( Track trk : recoTracks ) {
+ List<CalorimeterHit> hits = trkHitsMap.get(trk);
+ if(_debug>0) {
+ System.out.println("Track "+i+": #match hits="
+ +( hits!=null ? hits.size() : 0 ) );
+ }
+ trackHitRels.add( new TrackHitRelation(trk, hits) );
+ ++i;
+ }
+
+ // stores track-hit associations into event
+ event.put("TrkHitsSwimmer", trackHitRels, TrackHitRelation.class, 0);
+ }
+
+ private void reset() {
+ }
+
+ private void initialize(EventHeader event) {
+ CalHitMapMgr expert = CalHitMapMgr.getInstance();
+
+ // face of EM calorimeter
+ CylindricalCalorimeter embSubdet = (CylindricalCalorimeter)expert.getSubdetector(_embName);
+ _rhoMinEM = embSubdet.getInnerRadius();
+ Layering layers = embSubdet.getLayering();
+ int nlayers = layers.getLayerCount();
+ _layersEMB = new double[nlayers];
+ for(int i=0; i<nlayers; ++i) {
+ _layersEMB[i] = layers.getDistanceToLayerSensorMid(i);
+ }
+
+ CylindricalCalorimeter emeSubdet = (CylindricalCalorimeter)expert.getSubdetector(_emeName);
+ layers = emeSubdet.getLayering();
+ nlayers = layers.getLayerCount();
+ _layersEME = new double[nlayers];
+ for(int i=0; i<nlayers; ++i) {
+ _layersEME[i] = layers.getDistanceToLayerSensorMid(i);
+ }
+
+ // face of HAD calorimeter
+ CylindricalCalorimeter hdbSubdet = (CylindricalCalorimeter)expert.getSubdetector(_hdbName);
+ _rhoMinHAD = hdbSubdet.getInnerRadius();
+ layers = hdbSubdet.getLayering();
+ nlayers = layers.getLayerCount();
+ _layersHDB = new double[nlayers];
+ for(int i=0; i<nlayers; ++i) {
+ _layersHDB[i] = layers.getDistanceToLayerSensorMid(i);
+ }
+
+ CylindricalCalorimeter hdeSubdet = (CylindricalCalorimeter)expert.getSubdetector(_hdeName);
+ layers = hdeSubdet.getLayering();
+ nlayers = layers.getLayerCount();
+ _layersHDE = new double[nlayers];
+ for(int i=0; i<nlayers; ++i) {
+ _layersHDE[i] = layers.getDistanceToLayerSensorMid(i);
+ }
+
+ // setup a swimmer
+ double[] pos = {0,0,0};
+ double[] field = event.getDetector().getFieldMap().getField(pos);
+ _swimmer = new HelixSwimmer( field[2] );
+
+ // setup track-hit matchers
+ _embMatcher = new TrackHitMatcher(_embName, _layersEMB, 0, 0, _debug);
+ _emeMatcher = new TrackHitMatcher(_emeName, _layersEME, 0, 0, _debug);
+ _hdbMatcher = new TrackHitMatcher(_hdbName, _layersHDB, 0, 0, _debug);
+ _hdeMatcher = new TrackHitMatcher(_hdeName, _layersHDE, 0, 0, _debug);
+
+ _init = true;
+ }
+
+ //***** FIELDS ****
+
+ private boolean _init = false;
+ private double _rhoMinEM;
+ private double _rhoMinHAD;
+
+ private String _embName = "EcalBarrHits";
+ private String _hdbName = "HcalBarrHits";
+ private String _emeName = "EcalEndcapHits";
+ private String _hdeName = "HcalEndcapHits";
+ private double[] _layersEMB;
+ private double[] _layersEME;
+ private double[] _layersHDB;
+ private double[] _layersHDE;
+ private TrackHitMatcher _embMatcher;
+ private TrackHitMatcher _emeMatcher;
+ private TrackHitMatcher _hdbMatcher;
+ private TrackHitMatcher _hdeMatcher;
+
+ private HelixSwimmer _swimmer;
+ private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+}
CVSspam 0.2.8