lcsim/src/org/lcsim/recon/cluster/directedtree
diff -u -r1.3 -r1.4
--- DTreeAnalysis.java 20 Feb 2006 22:47:12 -0000 1.3
+++ DTreeAnalysis.java 15 Mar 2006 12:27:40 -0000 1.4
@@ -21,6 +21,7 @@
import org.lcsim.geometry.util.BaseIDDecoder;
import org.lcsim.event.Track;
import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.recon.ztracking.cheater.CheatTrack;
public class DTreeAnalysis extends Driver {
@@ -37,7 +38,7 @@
AIDA aida = AIDA.defaultInstance();
_tree = aida.tree();
_tf = aida.analysisFactory().createTupleFactory(_tree);
- _tree.mkdir("DirectedTree");
+ _tree.mkdir("DirectedTree");
_tree.cd("DirectedTree");
bookDirectedTreeNtuple();
_tree.cd("..");
@@ -80,6 +81,7 @@
processCollection(hdbName,hdbhitmap);
processCollection(hdeName,hdehitmap);
+ fillTrackMatchTuple();
finalizeTuple();
}
else {
@@ -97,7 +99,7 @@
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,long cellid}","evtno"};
+ "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);
@@ -106,11 +108,11 @@
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"};
+ 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,long cellid}","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"};
+ 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,long cellid}","evtno"};
Class[] columnClassesMChd = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
_tupleMChd = _tf.create("MCHD","Gen Clusters in HD",columnNamesMChd,columnClassesMChd);
@@ -130,8 +132,6 @@
Class[] trkMatchClasses = {Integer.TYPE,ITuple.class};
_tupleTrkMatch = _tf.create("TrkMatch","Track matching to calorimeter hits",trkMatchNames,trkMatchClasses);
}
- else{
- }
}
public void fillDirectedTreeNtuple(List<BasicCluster> recon,
@@ -160,13 +160,11 @@
= 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);
+// System.out.println("reco tag: "+tag+", #hits="+cellVec.size());
+ if(cellVec.size()==0) System.out.println("*** DTreeAna: 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);
+// System.out.println("putting <"+MyTools.printID(ihit.getCellID())+" -> "+tag+">");
+ recTagMap.put(ihit,tag);
double[] pos = ihit.getPosition();
double cdens = _loader.getDensity(ihit);
long jid = ihit.getCellID();
@@ -175,14 +173,6 @@
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;
@@ -194,10 +184,8 @@
cellFolder.fill(4,(float)cdens);
cellFolder.fill(5,tag);
cellFolder.fill(6,ily);
- if(_evtnum==223) System.out.println(_calType+", ily="+ily);
cellFolder.fill(7,iu);
cellFolder.fill(8,iv);
- cellFolder.fill(9,jid);
cellFolder.addRow();
}
}
@@ -217,21 +205,21 @@
}
if(_calType=="EM") _nClusEM = tag;
if(_calType=="HD") _nClusHD = tag;
-// System.out.println("calType="+_calType+", tag="+tag+", nrecoEM="+_nClusEM+", nrecoHD="+_nClusHD);
+// System.out.println("DTreeAna: 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 ) {
+
// loop over all hits on this MCcluster
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+">");
+ else System.out.println("***** DTreeAna: Invalid tag, cellid="+MyTools.printID(ihit.getCellID()));
double[] pos = ihit.getPosition();
double cdens = _loader.getDensity(ihit);
if(cdens==0) continue; // skip hits below energy cut
@@ -241,14 +229,6 @@
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]);
@@ -273,7 +253,6 @@
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);
@@ -283,132 +262,54 @@
mcCellFolder.fill(11,(float) phi);
mcCellFolder.fill(12,(float) theta);
mcCellFolder.fill(13,(float) timeStamp);
+ mcCellFolder.fill(14, ihit.getCellID() );
mcCellFolder.addRow();
-// System.out.println("EM Clus cells tag="+tag+", rtag="+rtag
-// +", id="+MyTools.printID(jid)
+// System.out.println("EM Clus cells tag="+tag+", rtag="+rtag
+// +", id="+MyTools.printID(jid)
// // +", pos="+pos[0]+" "+pos[1]+" "+pos[2]
-// +", dens="+cdens);
+// +", 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) {
+ ITuple mcClusFolder = null;
+ if(_calType=="EM") mcClusFolder = emmcClusFolder;
+ if(_calType=="HD") mcClusFolder = hdmcClusFolder;
+ if(mcClusFolder!=null) {
mcClusFolder.fill(0,(float) clust.getRawEnergy());
+ double[] pos = clust.getPosition();
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());
+ if(mcpart==null) {
+ // no MC contribution for this hit (e.g. from DigiSim)
+ mcClusFolder.fill(5, 999999);
+ }
+ else {
+ float charge = (float)mcpart.getCharge();
+ int particleID = mcpart.getType().getPDGID();
+ float mcE = (float)mcpart.getEnergy();
+ Hep3Vector mcMom = mcpart.getMomentum();
+ Hep3Vector vert = mcpart.getOrigin();
+ mcClusFolder.fill(5, particleID);
+ mcClusFolder.fill(6, charge);
+ mcClusFolder.fill(7, mcE);
+ mcClusFolder.fill(8, (float)mcMom.x());
+ mcClusFolder.fill(9,(float)mcMom.y());
+ mcClusFolder.fill(10,(float)mcMom.z());
+ mcClusFolder.fill(11,(float) vert.x());
+ mcClusFolder.fill(12,(float) vert.y());
+ mcClusFolder.fill(13,(float) vert.z());
+ }
mcClusFolder.addRow();
- }
- tag++;
}
+ tag++;
}
if(_calType=="EM") _nMCClusEM = tag;
if(_calType=="HD") _nMCClusHD = tag;
// System.out.println("calType="+_calType+", tag="+tag+", nMCEM="+_nMCClusEM+", nMCHD="+_nMCClusHD);
-
- // Fill tuples with track match info (swimmer)
- List<TrackHitRelation> matches = _evt.get(TrackHitRelation.class, "TrkHitsSwimmer");
- assert matches!=null
- : "DTreeAnalysis: no track-hit matching information available!";
-
- _tupleTrkMatch.fill(0, matches.size() );
-
- ITuple tracksTuple = _tupleTrkMatch.getTuple(1);
- for( TrackHitRelation tkhit : matches ) {
- Track trk = tkhit.getTrack();
- Hep3Vector pvec = new BasicHep3Vector( trk.getMomentum() );
- List<CalorimeterHit> hits = tkhit.getHits();
- if( hits == null ) continue;
-
- if(_debug>0) {
- System.out.println("Swimmer: Track |p|="+pvec.magnitude()
- +", #match hits="+(hits!=null ? hits.size() : 0));
- }
-
- MCParticle mcp = (MCParticle) ((ReconTrack)trk).getMCParticle();
- assert mcp!=null : "DTreeAnalysis: MCParticle info not available in <TrkHitsSwimmer> collection";
-
- float charge = (float)mcp.getCharge();
- int particleID = mcp.getType().getPDGID();
- float mcE = (float)mcp.getEnergy();
- Hep3Vector mcMom = mcp.getMomentum();
- Hep3Vector vert = mcp.getOrigin();
-
- tracksTuple.fill(0, particleID);
- tracksTuple.fill(1, charge);
- tracksTuple.fill(2, mcE);
- tracksTuple.fill(3, (float)mcMom.x());
- tracksTuple.fill(4, (float)mcMom.y());
- tracksTuple.fill(5, (float)mcMom.z());
- tracksTuple.fill(6, (float) vert.x());
- tracksTuple.fill(7, (float) vert.y());
- tracksTuple.fill(8, (float) vert.z());
-
- tracksTuple.fill(9, hits.size());
- ITuple hitsTuple = tracksTuple.getTuple(10);
- for(CalorimeterHit hit : hits) {
- hitsTuple.fill(0,(float) hit.getRawEnergy());
- double[] pos = hit.getPosition();
- hitsTuple.fill(1,(float) pos[0]);
- hitsTuple.fill(2,(float) pos[1]);
- hitsTuple.fill(3,(float) pos[2]);
- hitsTuple.fill(4,(float) _loader.getDensity(hit.getCellID()));
-
- decoder.setID( hit.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
- hitsTuple.fill(5, ily);
- hitsTuple.fill(6, iu);
- hitsTuple.fill(7, iv);
-
- double r = Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
- hitsTuple.fill(8, (float)r);
- double phi = Math.atan2(pos[1], pos[0]);
- if(phi<0) phi += 2 * Math.PI;
- hitsTuple.fill(9, (float)phi);
- double theta = Math.atan2(r, pos[2]);
- hitsTuple.fill(10, (float)theta);
-
- double timeStamp = 99999.;
- // loop over contributions to this hit
- SimCalorimeterHit iihit = (SimCalorimeterHit)hit;
- for(int j=0;j<iihit.getMCParticleCount();j++) {
- double tmpTime = (double) iihit.getContributedTime(j);
- if(tmpTime<timeStamp){
- timeStamp = tmpTime;
- }
- }
- hitsTuple.fill(11, (float)timeStamp);
- hitsTuple.fill(12, hit.getCellID());
- hitsTuple.addRow();
- }
- tracksTuple.addRow();
- }
- _tupleTrkMatch.addRow();
-
recTagMap.clear();
}
@@ -428,6 +329,7 @@
_tupleMChd.fill(3,_evtnum);
_tupleMChd.addRow();
// event
+ _tupleEvt.fill(0,_evtnum);
_tupleEvt.addRow();
}
@@ -526,6 +428,90 @@
return null;
}
+ private void fillTrackMatchTuple() {
+ // Fill tuples with track match info (swimmer)
+ List<TrackHitRelation> matches = _evt.get(TrackHitRelation.class, "TrkHitsSwimmer");
+ assert matches!=null
+ : "DTreeAnalysis: no track-hit matching information available!";
+
+ int nmatches = 0;
+ ITuple tracksTuple = _tupleTrkMatch.getTuple(1);
+ for( TrackHitRelation tkhit : matches ) {
+ Track trk = tkhit.getTrack();
+ Hep3Vector pvec = new BasicHep3Vector( trk.getMomentum() );
+ List<CalorimeterHit> hits = tkhit.getHits();
+ if( hits == null ) continue;
+
+ ++nmatches;
+ if(_debug>0) {
+ System.out.println("Swimmer: Track |p|="+pvec.magnitude()
+ +", #match hits="+(hits!=null ? hits.size() : 0));
+ }
+
+ MCParticle mcp = (MCParticle) ((CheatTrack)trk).getMCParticle();
+ assert mcp!=null : "DTreeAnalysis: MCParticle info not available in <TrkHitsSwimmer> collection";
+
+ float charge = (float)mcp.getCharge();
+ int particleID = mcp.getType().getPDGID();
+ float mcE = (float)mcp.getEnergy();
+ Hep3Vector mcMom = mcp.getMomentum();
+ Hep3Vector vert = mcp.getOrigin();
+
+ tracksTuple.fill(0, particleID);
+ tracksTuple.fill(1, charge);
+ tracksTuple.fill(2, mcE);
+ tracksTuple.fill(3, (float)mcMom.x());
+ tracksTuple.fill(4, (float)mcMom.y());
+ tracksTuple.fill(5, (float)mcMom.z());
+ tracksTuple.fill(6, (float) vert.x());
+ tracksTuple.fill(7, (float) vert.y());
+ tracksTuple.fill(8, (float) vert.z());
+
+ tracksTuple.fill(9, hits.size());
+
+ ITuple hitsTuple = tracksTuple.getTuple(10);
+ for(CalorimeterHit hit : hits) {
+ BaseIDDecoder dec = (BaseIDDecoder)hit.getIDDecoder();
+ hitsTuple.fill(0,(float) hit.getRawEnergy());
+ double[] pos = hit.getPosition();
+ hitsTuple.fill(1,(float) pos[0]);
+ hitsTuple.fill(2,(float) pos[1]);
+ hitsTuple.fill(3,(float) pos[2]);
+ hitsTuple.fill(4,(float) _loader.getDensity(hit.getCellID()));
+ dec.setID( hit.getCellID() );
+ int ily = dec.getValue(0); // layer
+ int iu = dec.getValue(3); // either theta, phi or x
+ int iv = dec.getValue(4); // either phi, z or y
+ hitsTuple.fill(5, ily);
+ hitsTuple.fill(6, iu);
+ hitsTuple.fill(7, iv);
+
+ double r = Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
+ hitsTuple.fill(8, (float)r);
+ double phi = Math.atan2(pos[1], pos[0]);
+ if(phi<0) phi += 2 * Math.PI;
+ hitsTuple.fill(9, (float)phi);
+ double theta = Math.atan2(r, pos[2]);
+ hitsTuple.fill(10, (float)theta);
+
+ double timeStamp = 99999.;
+ // loop over contributions to this hit
+ SimCalorimeterHit iihit = (SimCalorimeterHit)hit;
+ for(int j=0;j<iihit.getMCParticleCount();j++) {
+ double tmpTime = (double) iihit.getContributedTime(j);
+ if(tmpTime<timeStamp){
+ timeStamp = tmpTime;
+ }
+ }
+ hitsTuple.fill(11, (float)timeStamp);
+ hitsTuple.fill(12, hit.getCellID());
+ hitsTuple.addRow();
+ }
+ tracksTuple.addRow();
+ }
+ _tupleTrkMatch.fill(0, nmatches );
+ _tupleTrkMatch.addRow();
+ }
// *** FIELDS ***