lcsim/src/org/lcsim/recon/cluster/directedtree
diff -u -r1.2 -r1.3
--- DTreeAnalysis.java 9 Feb 2006 21:17:59 -0000 1.2
+++ DTreeAnalysis.java 20 Feb 2006 22:47:12 -0000 1.3
@@ -1,6 +1,7 @@
package org.lcsim.recon.cluster.directedtree;
import java.io.IOException;
+import java.text.NumberFormat;
import java.util.*;
import hep.aida.*;
import hep.physics.vec.Hep3Vector;
@@ -18,15 +19,19 @@
import org.lcsim.recon.cluster.util.BasicCluster;
import org.lcsim.recon.cluster.util.SortClustersBySize;
import org.lcsim.geometry.util.BaseIDDecoder;
+import org.lcsim.event.Track;
+import org.lcsim.mc.fast.tracking.ReconTrack;
public class DTreeAnalysis extends Driver {
private int _debug=0;
+ // constructor
public DTreeAnalysis() {
_runPar = RunControlParameters.getInstance();
_emEcut = _runPar.getEMmip() * _runPar.getEMthresh() / _runPar.getEMweight();
_hdEcut = _runPar.getHDmip() * _runPar.getHDthresh() / _runPar.getHDweight();
+ _format.setMinimumFractionDigits(12);
// aida stuff
AIDA aida = AIDA.defaultInstance();
@@ -40,6 +45,7 @@
public void process(EventHeader event) {
+ _evt = event;
_evtnum = event.getEventNumber();
if(_debug>0) {
System.out.println("Starting DTreeAnalysis, event # "+_evtnum);
@@ -66,13 +72,13 @@
if(_runPar.ClusterSeparately()){
_calType = "EM";
_nClusEM = 0; _nMCClusEM = 0;
- processCollection(event,embName,embhitmap);
- processCollection(event,emeName,emehitmap);
+ processCollection(embName,embhitmap);
+ processCollection(emeName,emehitmap);
_calType = "HD";
_nClusHD = 0; _nMCClusHD = 0;
- processCollection(event,hdbName,hdbhitmap);
- processCollection(event,hdeName,hdehitmap);
+ processCollection(hdbName,hdbhitmap);
+ processCollection(hdeName,hdehitmap);
finalizeTuple();
}
@@ -91,12 +97,12 @@
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"};
+ "emCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTag,int ily,int iz,int iphi,long cellid}","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"};
+ "hdCellFolder={float cellE,float cellX,float cellY,float cellZ,float cellD,int cellTag,int ily,int iz,int iphi,long cellid}","evtno"};
Class[] columnClassesHD = {Integer.TYPE,ITuple.class,ITuple.class,Integer.TYPE};
_tupleHD = _tf.create("HD","Recon Clusters",columnNamesHD,columnClassesHD);
@@ -107,6 +113,22 @@
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);
+
+ /* commented out for now... all hits are already in
+ * MCXXcellFolder, unless a minimum cluster size is
+ * required
+// // add tuple with *all* hits above threshold
+// String[] columnNamesHits = {"nhits",
+// "hitFolder={float hitE,float hitX,float hitY,float hitZ,float hitD,int hitTag,int ily,int iz,int iphi,float r,float phi,float theta,float time}","evtno"};
+// Class[] columnClassesHits = {Integer.TYPE,ITuple.class,Integer.TYPE};
+// _tupleHits = _tf.create("Hits","Cal hits in event",
+// columnNamesHits,columnClassesHits);
+ */
+
+ // add tuple for track-hit matching
+ String[] trkMatchNames = {"ntrkmatch","trkMatches={int pid,float chrg,float e,float px,float py,float pz,float vx,float vy,float vz,int nhits,ITuple hits={float cellE,float cellX,float cellY,float cellZ,float cellD,int ily,int iz,int iphi,float r,float phi,float theta,float time,long cellid}}"};
+ Class[] trkMatchClasses = {Integer.TYPE,ITuple.class};
+ _tupleTrkMatch = _tf.create("TrkMatch","Track matching to calorimeter hits",trkMatchNames,trkMatchClasses);
}
else{
}
@@ -175,6 +197,7 @@
if(_evtnum==223) System.out.println(_calType+", ily="+ily);
cellFolder.fill(7,iu);
cellFolder.fill(8,iv);
+ cellFolder.fill(9,jid);
cellFolder.addRow();
}
}
@@ -304,6 +327,88 @@
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();
}
@@ -323,11 +428,10 @@
_tupleMChd.fill(3,_evtnum);
_tupleMChd.addRow();
// event
- _tupleEvt.fill(0,_evtnum);
_tupleEvt.addRow();
}
- private void processCollection(EventHeader event, String colName,
+ private void processCollection(String colName,
Map<Long,CalorimeterHit> hitmap)
{
if(hitmap.size()==0) return;
@@ -336,7 +440,7 @@
// get NearestNeighbor clusters
List<BasicCluster> nnClusColl = new ArrayList<BasicCluster>();
try {
- nnClusColl = event.get(BasicCluster.class, colName+"NNClusters");
+ nnClusColl = _evt.get(BasicCluster.class, colName+"NNClusters");
}
catch(IllegalArgumentException x) {
// ignore
@@ -345,7 +449,7 @@
// get reco clusters
List<BasicCluster> recoClusColl = new ArrayList<BasicCluster>();
try {
- recoClusColl = event.get(BasicCluster.class, colName+"DTreeClusters");
+ recoClusColl = _evt.get(BasicCluster.class, colName+"DTreeClusters");
}
catch(IllegalArgumentException x) {
// ignore
@@ -353,7 +457,7 @@
// get MC clusters
List<BasicCluster> mcClusColl
- = event.get(BasicCluster.class, colName+"CheatClusters");
+ = _evt.get(BasicCluster.class, colName+"CheatClusters");
Collections.sort( mcClusColl, new SortClustersBySize() );
if(_debug>0) {
@@ -422,6 +526,10 @@
return null;
}
+
+ // *** FIELDS ***
+
+ private static EventHeader _evt;
private static Map<String,Readout> _roMap;
private static RunControlParameters _runPar;
private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
@@ -444,9 +552,11 @@
private ITuple hdClusFolder,hdCellFolder;
private ITuple emmcClusFolder,emmcCellFolder;
private ITuple hdmcClusFolder,hdmcCellFolder;
+ private ITuple _tupleTrkMatch;
private ClusterBuilder clusBuilder = new ClusterBuilder();
private int layerIndex, uIndex, vIndex;
private double _emEcut, _hdEcut;
private int _nClusEM,_nMCClusEM,_nClusHD,_nMCClusHD;
+ private NumberFormat _format = NumberFormat.getInstance();
}
