lcsim/src/org/lcsim/recon/cluster/density
diff -N TrackMatchingDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackMatchingDriver.java 27 Dec 2005 16:08:36 -0000 1.1
@@ -0,0 +1,237 @@
+package org.lcsim.recon.cluster.density;
+
+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 2005/12/27 16:08:36 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();
+}
lcsim/src/org/lcsim/recon/cluster/density
diff -N TrackMatching.java
--- TrackMatching.java 12 Dec 2005 05:32:27 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,635 +0,0 @@
-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.segmentation.BarrelCylinderSegmentationBase;
-import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
-import org.lcsim.geometry.layer.Layering;
-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;
-
-/**
- * A driver for track matching using the swimmer.
- *
- * @author Guilherme Lima
- * @version $Id: TrackMatching.java,v 1.2 2005/12/12 05:32:27 lima Exp $
- */
-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;
-
- int nhitsInLayerEmb = _loader.getHitCountInLayer(_embName,i);
- int nhitsInLayerEme = _loader.getHitCountInLayer(_emeName,i);
- int nhitsInLayer = nhitsInLayerEmb + nhitsInLayerEme;
- System.out.println("#hits in layer: "+nhitsInLayerEmb+" "+nhitsInLayerEme+" "+nhitsInLayer);
- if(nhitsInLayer<1) continue;
-
- 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
- LoadMyCalorimeterHit aux = LoadMyCalorimeterHit.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 +/-2 neighbors
- long[] neighs2 = _segm.getNeighbourIDs(0,2,2);
- nmatch = 0;
- nbad = 0;
- if(perfmatch) {
- if( isContributor(ipart,ihit) ) ++nmatch;
- else ++nbad;
- }
- for(int j=0; j<neighs2.length; ++j) {
- CalorimeterHit jhit = embhitmap.get( neighs2[j] );
- if( jhit!=null ) {
- if( isContributor(ipart,jhit) ) ++nmatch;
- else ++nbad;
- }
- }
-// System.out.println("good022="+nmatch+", bad022="+nbad);
- _aida.cloud1D("trkMatch 022").fill(nmatch);
- if(nbad>0) _aida.cloud1D("badMatch 022").fill(nbad);
-
- // 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);
- }
- }
-
-// 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) {
- LoadMyCalorimeterHit expert = LoadMyCalorimeterHit.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 HAD 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;
-}