lcsim/src/org/lcsim/recon/cluster/density
diff -N TrackMatching.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackMatching.java 7 Dec 2005 18:57:20 -0000 1.1
@@ -0,0 +1,624 @@
+package org.lcsim.recon.cluster.density;
+
+import java.io.IOException;
+import java.util.*;
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.geometry.util.CalorimeterIDDecoder;
+import org.lcsim.geometry.compact.Subdetector;
+// import org.lcsim.geometry.compact.Readout;
+import org.lcsim.geometry.segmentation.BarrelCylinderSegmentationBase;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.layer.Layering;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.SortClustersBySize;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.util.swim.HelixSwim;
+
+public class TrackMatching extends Driver {
+
+ int _debug = 1; // debug level, 0 for no printout
+ boolean _fillTuple = false;
+
+ public TrackMatching() throws IOException {
+// _runPar = RunControlParameters.getInstance();
+ _aida = AIDA.defaultInstance();
+ _tree = _aida.tree();
+ _tf = _aida.analysisFactory().createTupleFactory(_tree);
+ _tree.mkdir("TrkClusMatch");
+ _tree.cd("TrkClusMatch");
+ if(_fillTuple) bookTrkClusMatchNtuple();
+ _tree.cd("..");
+ }
+
+ public void process(EventHeader event) {
+
+ System.out.println("******** Track to Cluster matching *******");
+ _tree.cd("TrkClusMatch");
+
+ if(!_init) initialize(event);
+ this.reset();
+
+ // Load hit collections
+ embhitmap = _loader.getCollHitMap(_embName);
+ hdbhitmap = _loader.getCollHitMap(_hdbName);
+// emehitmap = _loader.getCollHitMap(_emeName);
+// hdehitmap = _loader.getCollHitMap(_hdeName);
+
+ _evtnum = event.getEventNumber();
+
+ List<MCParticle> evtParticles = event.getMCParticles();
+ for( MCParticle ipart : evtParticles ) {
+ int id = ipart.getPDGID();
+ int absid = Math.abs(id);
+ double energy = ipart.getEnergy();
+ Hep3Vector orig = ipart.getOrigin();
+ Hep3Vector endp = ipart.getEndPoint();
+ double rhoOrig = Math.sqrt(orig.x()*orig.x() + orig.y()*orig.y());
+ double rhoEndp = Math.sqrt(endp.x()*endp.x() + endp.y()*endp.y());
+
+ // save particles which reach EM calorimeter face
+ System.out.println("particle: "+id+" "+ipart.getType().getName()
+ +", discard="+this.discardParticleType(ipart) );
+ if( this.discardParticleType(ipart) ) continue;
+ boolean origBeforeEM = rhoOrig<_rhoMinEM && Math.abs(orig.z())<_zMinEM;
+ boolean endpBeforeEM = rhoEndp<_rhoMinEM && Math.abs(endp.z())<_zMinEM;
+ System.out.println("particle: "+id+" "+ipart.getType().getName()
+ +", E="+ipart.getEnergy()
+ +", genStat="+ipart.getGeneratorStatus()
+ +", rho="+rhoOrig+" "+rhoEndp
+ +", z="+orig.z()+" "+endp.z()
+ +", origBeforeEM="+origBeforeEM
+ +", endpBeforeEM="+endpBeforeEM
+ );
+
+ if ( origBeforeEM && !endpBeforeEM ) {
+
+ if (ipart.getGeneratorStatus() == 0)
+ _aida.cloud1D("rhoOrig, genStat=0").fill(rhoOrig);
+
+ if(id==22) _photons.add(ipart);
+ else if(ipart.getCharge()==0) _neutrals.add(ipart);
+ else {
+ System.out.println("Charged particle saved!");
+ if(ipart.getEnergy()>0.5) _chgParts.add(ipart);
+ }
+ }
+ }
+
+ System.out.println("# parts: photons="+_photons.size()
+ +", neutrals="+_neutrals.size()
+ +", chgTracks="+_chgParts.size() );
+
+ // loop over particles matching them to clusters
+ for(MCParticle ipart : _chgParts) {
+
+ System.out.println("*** New particle: "+ipart.getType().getName()
+ +", E="+ipart.getEnergy());
+
+ // Swim the particle to each sensitive layer
+ boolean curling = false;
+ int icharge = (int)ipart.getCharge();
+ _swimmer.setTrack( ipart.getMomentum(), ipart.getOrigin(), icharge );
+ for(int i=0; i<_layersEMB.length; ++i) {
+ if(curling) 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) {
+ curling = true;
+ System.out.println("Track seems to be curling. Break.");
+ break;
+ }
+ }
+ System.out.println("Swimmer: layer="+i+", pos="+pos);
+
+ // find cell containing track-cylinder intersection
+ CalHitMapMgr aux = CalHitMapMgr.getInstance();
+ Subdetector _calsub = aux.getSubdetector( "EcalBarrHits" );
+ BarrelCylinderSegmentationBase _segm = (BarrelCylinderSegmentationBase)_calsub.getReadout().getSegmentation();
+
+ long cellid = _segm.findCellContainingXYZ( pos );
+ CalorimeterHit ihit = embhitmap.get(cellid);
+// System.out.println("from findXYZ: cellid="+MyTools.printID(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="+embhitmap.size());
+
+ boolean perfmatch = false;
+ SimCalorimeterHit mchit = null;
+ if(ihit!=null) {
+ // perfect track match
+ perfmatch = true;
+ _aida.cloud1D("NN energy").fill(((SimCalorimeterHit)ihit).getRawEnergy());
+ _aida.cloud1D("trkMatch 000").fill(1);
+ }
+ else {
+ _aida.cloud1D("trkMatch 000").fill(0);
+ }
+
+ // look for hits in central cells and neighbors
+ long[] neighs1 = _segm.getNeighbourIDs(0,1,1);
+ // use different counters for good and bad match
+ int nmatch = 0;
+ int nbad = 0;
+ if( perfmatch ) {
+ if( isContributor(ipart,ihit) ) ++nmatch;
+ else ++nbad;
+ }
+ for(int j=0; j<neighs1.length; ++j) {
+ CalorimeterHit jhit = embhitmap.get( neighs1[j] );
+ if(jhit!=null) {
+ _aida.cloud1D("NN energy").fill(((SimCalorimeterHit)jhit).getRawEnergy());
+ if( isContributor(ipart,jhit) ) ++nmatch;
+ else ++nbad;
+ }
+ }
+// System.out.println("good011="+nmatch+", bad011="+nbad);
+ _aida.cloud1D("trkMatch 011").fill(nmatch);
+ if(nbad>0) {
+ _aida.cloud1D("badMatch 011").fill(nbad);
+ System.out.println("a bad match: event# "+_evtnum);
+ }
+
+ // look for +/-3 neighbors
+ long[] neighs3 = _segm.getNeighbourIDs(0,3,3);
+ nmatch = 0;
+ nbad = 0;
+ if(perfmatch) {
+ if( isContributor(ipart,ihit) ) ++nmatch;
+ else ++nbad;
+ }
+ for(int j=0; j<neighs3.length; ++j) {
+ CalorimeterHit jhit = embhitmap.get( neighs3[j] );
+ if( jhit!=null ) {
+ if( isContributor(ipart,jhit) ) ++nmatch;
+ else ++nbad;
+ }
+ }
+// System.out.println("good033="+nmatch+", bad033="+nbad);
+ _aida.cloud1D("trkMatch 033").fill(nmatch);
+ if(nbad>0) _aida.cloud1D("badMatch 033").fill(nbad);
+
+ // look for +/-5 neighbors
+ long[] neighs5 = _segm.getNeighbourIDs(0,5,5);
+ nmatch = 0;
+ nbad = 0;
+ if(perfmatch) {
+ if( isContributor(ipart,ihit) ) ++nmatch;
+ else ++nbad;
+ }
+ for(int j=0; j<neighs5.length; ++j) {
+ CalorimeterHit jhit = embhitmap.get( neighs5[j] );
+ if( jhit!=null ) {
+ if( isContributor(ipart,jhit) ) ++nmatch;
+ else ++nbad;
+ }
+ }
+// System.out.println("good055="+nmatch+", bad055="+nbad);
+ _aida.cloud1D("trkMatch 055").fill(nmatch);
+ if(nbad>0) _aida.cloud1D("badMatch 055").fill(nbad);
+ }
+ }
+
+// if(_runPar.ClusterSeparately()){
+// processCollection(event,_embName);
+// // processCollection(event,_emeName);
+
+// processCollection(event,_hdbName);
+// // processCollection(event,_hdeName);
+// }
+// else {
+// assert false : "Sorry, single-pass clustering unavailable for now.";
+// }
+
+ _tree.cd("..");
+ }
+
+ public boolean discardParticleType(MCParticle ipart) {
+ int id = ipart.getPDGID();
+ int absid = Math.abs(id);
+ // discard quarks
+ if( absid<10 ) return true;
+ // discard neutrinos
+ if(absid==12 || absid==14 || absid==16) return true;
+ // discard strings
+ if(id>=90 && id<=100) return true;
+ return false;
+ }
+
+ public void bookTrkClusMatchNtuple(){
+// 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 fillTrkClusMatchNtuple(List<BasicCluster> recon,
+ List<BasicCluster> mctruth)
+ {
+// if(_calType=="EM"){
+// emClusFolder = _tupleEM.getTuple(1);
+// emCellFolder = _tupleEM.getTuple(2);
+// emmcClusFolder = _tupleMCem.getTuple(1);
+// emmcCellFolder = _tupleMCem.getTuple(2);
+// }
+
+// if(_calType=="HD"){
+// hdClusFolder = _tupleHD.getTuple(1);
+// hdCellFolder = _tupleHD.getTuple(2);
+// hdmcClusFolder = _tupleMChd.getTuple(1);
+// hdmcCellFolder = _tupleMChd.getTuple(2);
+// }
+
+// int tag = 0;
+// Map<CalorimeterHit,Integer> recTagMap
+// = new HashMap<CalorimeterHit,Integer>();
+// for(BasicCluster clust : recon ) {
+// List<CalorimeterHit> cellVec = clust.getCalorimeterHits();
+// // System.out.println(tag+" "+cellVec.size());
+// if(cellVec.size()==0.) System.out.println("Damn 3");
+// // for(int i=0;i<cellVec.size();i++){
+// // MyCalorimeterHit myhit = (MyCalorimeterHit) cellVec.get(i);
+// for(CalorimeterHit ihit : cellVec) {
+// Integer tagObject = new Integer(tag);
+// // System.out.println("putting <"+myhit+"->"+tagObject+">");
+// recTagMap.put(ihit,tagObject);
+// double[] pos = ihit.getPosition();
+// double cdens = _loader.getDensity(ihit);
+// if(_calType=="EM"){
+// emCellFolder.fill(0,(float) ihit.getRawEnergy());
+// emCellFolder.fill(1,(float)pos[0]);
+// emCellFolder.fill(2,(float)pos[1]);
+// emCellFolder.fill(3,(float)pos[2]);
+// emCellFolder.fill(4,(float) cdens);
+// emCellFolder.fill(5,tag);
+// emCellFolder.addRow();
+// }
+// if(_calType=="HD"){
+// hdCellFolder.fill(0,(float) ihit.getRawEnergy());
+// hdCellFolder.fill(1,(float)pos[0]);
+// hdCellFolder.fill(2,(float)pos[1]);
+// hdCellFolder.fill(3,(float)pos[2]);
+// hdCellFolder.fill(4,(float) cdens);
+// hdCellFolder.fill(5,tag);
+// hdCellFolder.addRow();
+// }
+// }
+// ITuple clusFolder = null;
+// if(_calType=="EM") clusFolder = emClusFolder;
+// if(_calType=="HD") clusFolder = hdClusFolder;
+// if(clusFolder!=null) {
+// Hep3Vector pos = new BasicHep3Vector( clust.getPosition() );
+// clusFolder.fill(0,(float) clust.getEnergy());
+// clusFolder.fill(1,(float) pos.x());
+// clusFolder.fill(2,(float) pos.y());
+// clusFolder.fill(3,(float) pos.z());
+// clusFolder.fill(4, clust.getSize());
+// clusFolder.addRow();
+// }
+// // if(_calType=="HD"){
+// // hdClusFolder.fill(0,(float) clust.getEnergy());
+// // hdClusFolder.fill(1,(float) clust._xpos);
+// // hdClusFolder.fill(2,(float) clust._ypos);
+// // hdClusFolder.fill(3,(float) clust._zpos);
+// // hdClusFolder.fill(4, clust.getSize() );
+// // hdClusFolder.addRow();
+// // }
+// tag++;
+// }
+
+// tag = 0;
+// for( BasicCluster clust : mctruth ) {
+// for( CalorimeterHit ihit : clust.getCalorimeterHits() ) {
+// Integer tagObject = recTagMap.get(ihit);
+// // System.out.println("getting: "+ihit+" -> "+tagObject);
+
+// int rtag = -1;
+// if(tagObject!=null) rtag = tagObject.intValue();
+// double[] pos = ihit.getPosition();
+// double cdens = _loader.getDensity(ihit);
+// if(cdens==0) continue; // skip hits below energy cut
+// long jid = ihit.getCellID();
+// int ily = MyTools.getLayer(jid);
+// int iz = MyTools.getThetaBin(jid);
+// int iphi = MyTools.getPhiBin(jid);
+// if(_calType=="EM") {
+// emmcCellFolder.fill(0,(float) ihit.getRawEnergy());
+// emmcCellFolder.fill(1,(float)pos[0]);
+// emmcCellFolder.fill(2,(float)pos[1]);
+// emmcCellFolder.fill(3,(float)pos[2]);
+// emmcCellFolder.fill(4,(float) cdens);
+// emmcCellFolder.fill(5,tag);
+// emmcCellFolder.fill(6,rtag);
+// emmcCellFolder.fill(7,ily);
+// emmcCellFolder.fill(8,iz);
+// emmcCellFolder.fill(9,iphi);
+// emmcCellFolder.addRow();
+// // System.out.println("EM Clus cells tag="+tag+", rtag="+rtag
+// // +", id="+MyTools.printID(jid)
+// // // +", pos="+pos[0]+" "+pos[1]+" "+pos[2]
+// // +", dens="+cdens);
+// }
+// if(_calType=="HD") {
+// hdmcCellFolder.fill(0,(float) ihit.getRawEnergy());
+// hdmcCellFolder.fill(1,(float)pos[0]);
+// hdmcCellFolder.fill(2,(float)pos[1]);
+// hdmcCellFolder.fill(3,(float)pos[2]);
+// hdmcCellFolder.fill(4,(float) cdens);
+// hdmcCellFolder.fill(5,tag);
+// hdmcCellFolder.fill(6,rtag);
+// hdmcCellFolder.fill(7,ily);
+// hdmcCellFolder.fill(8,iz);
+// hdmcCellFolder.fill(9,iphi);
+// hdmcCellFolder.addRow();
+// // if(tag==7){
+// // System.out.println("HD Clus cells: id="+MyTools.printID(jid)
+// // +", pos="+pos[0]+" "+pos[1]+" "+pos[2]
+// // +", dens="+cdens);
+// // }
+// }
+// }
+// MCParticle mcpart = getMCParticleInCluster(clust);
+
+// float charge = (float)mcpart.getCharge();
+// int particleID = mcpart.getType().getPDGID();
+// float mcE = (float)mcpart.getEnergy();
+// Hep3Vector mcMom = mcpart.getMomentum();
+// double[] pos = clust.getPosition();
+// if(_calType=="EM"){
+// emmcClusFolder.fill(0,(float) clust.getEnergy());
+// emmcClusFolder.fill(1, (float)pos[0] );
+// emmcClusFolder.fill(2, (float)pos[1] );
+// emmcClusFolder.fill(3, (float)pos[2] );
+// emmcClusFolder.fill(4, clust.getSize() );
+// emmcClusFolder.fill(5, particleID);
+// emmcClusFolder.fill(6, charge);
+// emmcClusFolder.fill(7, mcE);
+// emmcClusFolder.fill(8, (float)mcMom.x());
+// emmcClusFolder.fill(9,(float)mcMom.y());
+// emmcClusFolder.fill(10,(float)mcMom.z());
+// emmcClusFolder.addRow();
+// }
+// if(_calType=="HD") {
+// hdmcClusFolder.fill(0, (float)clust.getEnergy() );
+// hdmcClusFolder.fill(1, (float)pos[0] );
+// hdmcClusFolder.fill(2, (float)pos[1] );
+// hdmcClusFolder.fill(3, (float)pos[2] );
+// hdmcClusFolder.fill(4, clust.getSize() );
+// hdmcClusFolder.fill(5, particleID);
+// hdmcClusFolder.fill(6, (int) charge);
+// hdmcClusFolder.fill(7, mcE);
+// hdmcClusFolder.fill(8, (float)mcMom.x());
+// hdmcClusFolder.fill(9, (float)mcMom.y());
+// hdmcClusFolder.fill(10, (float)mcMom.z());
+// hdmcClusFolder.addRow();
+// }
+// tag++;
+// }
+
+// if(_calType=="EM"){
+// _tupleEM.fill(0,recon.size());
+// _tupleEM.fill(3,_evtnum);
+// _tupleEM.addRow();
+// _tupleMCem.fill(0,mctruth.size());
+// _tupleMCem.fill(3,_evtnum);
+// _tupleMCem.addRow();
+// _tupleEvt.fill(0,_evtnum);
+// _tupleEvt.addRow();
+// }
+// if(_calType=="HD"){
+// _tupleHD.fill(0,recon.size());
+// _tupleHD.fill(3,_evtnum);
+// _tupleHD.addRow();
+// _tupleMChd.fill(0,mctruth.size());
+// _tupleMChd.fill(3,_evtnum);
+// _tupleMChd.addRow();
+// }
+
+// recTagMap.clear();
+ }
+
+ private void processCollection(EventHeader event, String colName) {
+
+ Map<Long,CalorimeterHit> hitmap = _loader.getCollHitMap(colName);
+
+ // get reco clusters
+ List<BasicCluster> recoClusColl
+ = event.get(BasicCluster.class, colName+"DTreeClusters");
+
+ // get MC clusters
+ List<BasicCluster> mcClusColl
+ = event.get(BasicCluster.class, colName+"Clusters");
+ Collections.sort( mcClusColl, new SortClustersBySize() );
+
+ System.out.println("*** TrackMatching: colName="+colName
+ +", # of MCclusters: "+mcClusColl.size());
+
+ // Analyze MC clusters
+ int mchits=0;
+ for( BasicCluster iclus : mcClusColl ) {
+ double pos[] = iclus.getPosition();
+ List<CalorimeterHit> iHits = iclus.getCalorimeterHits();
+ mchits += iHits.size(); // includes hits below threshold
+ }
+ int nclusCheat = mcClusColl.size();
+
+ System.out.println(" Comparing cheat clusters: "
+ +", cheater="+nclusCheat);
+ System.out.println(" Comparing total #hits: "
+ +" in event: "+hitmap.size()
+ +", in MCClusters="+mchits );
+
+ int nhits = hitmap.size();
+ _aida.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);
+
+ if(_fillTuple) fillTrkClusMatchNtuple( recoClusColl, mcClusColl );
+ }
+
+ private boolean isContributor(MCParticle part, CalorimeterHit hit) {
+ if(hit==null) return false;
+ SimCalorimeterHit mchit = (SimCalorimeterHit)hit;
+ if(mchit==null) return false;
+ int nmc = mchit.getMCParticleCount();
+ for(int imc=0; imc<nmc; ++imc) {
+ if(mchit.getMCParticle(imc) == part) return true;
+ }
+ return false;
+ }
+
+ 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 void reset() {
+ _photons.clear();
+ _neutrals.clear();
+ _chgParts.clear();
+ }
+
+ 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);
+ _zMinEM = emeSubdet.getZMin();
+ layers = emeSubdet.getLayering();
+ nlayers = layers.getLayerCount();
+ _layersEME = new double[nlayers];
+ for(int i=0; i<nlayers; ++i) {
+ _layersEME[i] = layers.getDistanceToLayerSensorMid(i);
+ }
+
+ // face of EM calorimeter
+ CylindricalCalorimeter hdbSubdet = (CylindricalCalorimeter)expert.getSubdetector(_hdbName);
+ 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]);
+
+ _init = true;
+ }
+
+ //***** FIELDS ****
+
+ private String _embName = "EcalBarrHits";
+ private String _hdbName = "HcalBarrHits";
+ private String _emeName = "EcalEndcapHits";
+ private String _hdeName = "HcalEndcapHits";
+
+// private static Map<String,Readout> _roMap;
+// private static RunControlParameters _runPar;
+ private LoadMyCalorimeterHit _loader = LoadMyCalorimeterHit.getInstance();
+// private Map<Long,Long> _ParentMap;
+ private Map<Long,CalorimeterHit> embhitmap,emehitmap,hdbhitmap,hdehitmap;
+// private String _distType,_calType;
+// private int _nLyr;
+// private int _nZ;
+// private int _nPhi;
+ private int _evtnum;
+
+ private ITree _tree;
+ private ITupleFactory _tf;
+// private ITuple _tupleEM;
+// private ITuple _tupleHD;
+// private ITuple _tupleMCem,_tupleMChd;
+// private ITuple _tupleEvt;
+// private ITuple emClusFolder,emCellFolder;
+// private ITuple hdClusFolder,hdCellFolder;
+// private ITuple emmcClusFolder,emmcCellFolder;
+// private ITuple hdmcClusFolder,hdmcCellFolder;
+// private ClusterBuilder clusBuilder = new ClusterBuilder();
+
+ private boolean _init = false;
+ private double _rhoMinEM, _zMinEM;
+
+ private List<MCParticle> _photons = new ArrayList<MCParticle>();
+ private List<MCParticle> _chgParts = new ArrayList<MCParticle>();
+ private List<MCParticle> _neutrals = new ArrayList<MCParticle>();
+ private HelixSwimmer _swimmer;
+
+ private double[] _layersEMB;
+ private double[] _layersEME;
+ private double[] _layersHDB;
+ private double[] _layersHDE;
+ private AIDA _aida;
+}