lcsim/src/org/lcsim/contrib/SteveMagill
diff -N LEPhoDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LEPhoDriver.java 20 Aug 2008 20:21:53 -0000 1.1
@@ -0,0 +1,200 @@
+package org.lcsim.contrib.SteveMagill;
+
+// Photon Finder
+// Input : EM Clusters (DT, Cone)
+// 1) test clusters for proximity to track - if close enough, eliminate from photon list
+// default is 0.02 in track-cluster distance (theta-phi)
+// 2) H-Matrix test on most remaining clusters
+// Requirements : a) minimum of 25 hit pixels in cluster
+// b) minimum of 500 MeV energy (calibrated)
+// c) 1st layer hit < 8 (~5 X0)
+// Re-cluster hits with NN 1111 to define cluster core
+// if 1 NN cluster is defined and min number cells is 25, then apply H-Matrix test
+// else use original cluster and apply H-Matrix test
+// Minimum Log Diagonal Chisquared Probability is -9.9
+// Passing clusters are photons
+// 3) for remaining clusters :
+// if 1st layer < 8, photon if min energy > 100 MeV (calibrated)
+// else photon if :
+// no hadron cluster within 0.02 (same as above track-cluster distance)
+// and min energy > 100 MeV (calibrated)
+// output : list of photon clusters arranged by subdetector, list of photon reconstructed particles
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.util.Driver;
+import org.lcsim.util.swim.*;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.ITree;
+import org.lcsim.geometry.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.recon.ztracking.cheater.*;
+import org.lcsim.recon.cluster.fixedcone.*;
+import org.lcsim.recon.cluster.analysis.*;
+import org.lcsim.recon.emid.hmatrix.HMatrix;
+import org.lcsim.recon.emid.hmatrix.HMatrixTask;
+import org.lcsim.recon.emid.hmatrix.HMatrixBuilder;
+import org.lcsim.recon.emid.hmatrix.HMatrixConditionsConverter;
+import org.lcsim.math.chisq.ChisqProb;
+import org.lcsim.recon.cluster.nn.NearestNeighborCluster;
+import org.lcsim.spacegeom.*;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.detector.DetectorIdentifierHelper;
+import org.lcsim.detector.identifier.*;
+
+public class LEPhoDriver extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ private int _mincells;
+ private double _minE;
+ private double _maxE;
+ private int _minlayer;
+ private int _nLayers;
+ private ITree _tree;
+ private boolean _initialized = false;
+ private boolean phoD = false;
+ private String _clusname;
+ private String _oclname;
+
+ public LEPhoDriver(int mincells, double minE, double maxE, int minlayer)
+ {
+ _mincells = mincells;
+ _minE = minE;
+ _maxE = maxE;
+ _minlayer = minlayer;
+ }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event); // executes all added drivers
+
+ if(!_initialized)
+ {
+ CylindricalCalorimeter calsubB = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMBarrel");
+ _nLayers = calsubB.getLayering().getLayerCount();
+// System.out.println("found "+_nLayers+" layers in the EMBarrel");
+ _initialized = true;
+ }
+
+ // loop over all EM clusters
+ List<BasicCluster> LEPhoclusters = new ArrayList<BasicCluster>(); // pho clusters
+ try
+ {
+ List<BasicCluster> DTclusters = event.get(BasicCluster.class,_clusname);
+
+ for (BasicCluster dtclus : DTclusters)
+ {
+ double ecp[] = dtclus.getPosition();
+ double ecpx = ecp[0];
+ double ecpy = ecp[1];
+ double ecpz = ecp[2];
+ double ecR = Math.sqrt(ecpx*ecpx+ecpy*ecpy);
+ double ebclth = Math.atan(ecR/ecpz);
+ if (ebclth<0) ebclth+=Math.PI;
+ double ebclph = Math.atan2(ecpy,ecpx);
+ if (ebclph<0) ebclph+=2*Math.PI;
+ if (phoD) aida.cloud1D("Theta of LEPho Cand Clusters").fill(ebclth);
+ if (phoD) aida.cloud1D("Phi of LEPho Cand Clusters").fill(ebclph);
+
+ // get first layer info for this cluster
+ double[] layE = layerEnergies(dtclus);
+ int layint = 0;
+ for(int j=0; j<layE.length; ++j)
+ {
+ layint = j;
+ if (layE[j] > 0.02) break;
+ }
+ if (phoD) aida.cloud1D("Layer of first LEPho cand interaction").fill(layint);
+
+ // Now, first check interaction layer - must be >0 and < 9 or probably not photon - get out
+ if (layint<1 || layint>_minlayer) continue;
+ // check maxE first
+ if (dtclus.getEnergy()>_maxE) continue;
+ // now check minE, maxE, minnhits
+ if (dtclus.getEnergy()>_minE || dtclus.getCalorimeterHits().size()>_mincells)
+ {
+ LEPhoclusters.add(dtclus);
+ } // end of minE, mincell requirement loop
+ } // end of cluster loop
+ }
+ catch (java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No clusters in LEPho finder");
+ }
+ event.put(_oclname,LEPhoclusters);
+
+ // add both types of clusters here
+// List<BasicCluster> phoclusters = new ArrayList<BasicCluster>();
+// List<BasicCluster> hmphos = event.get(BasicCluster.class,"HMPhotonClusters");
+// for (BasicCluster phcl : LEPhoclusters)
+// {
+// phoclusters.add(phcl);
+// }
+// for (BasicCluster hmpho : hmphos)
+// {
+// phoclusters.add(hmpho);
+// }
+// event.put("FPhotonClusters",phoclusters);
+
+// if (phoD)
+// {
+// aida.cloud1D("Number of Photon Clusters").fill(phoclusters.size());
+// double PFAPhoE = 0;
+// for (BasicCluster pho : phoclusters)
+// {
+// PFAPhoE += pho.getEnergy();
+// aida.cloud1D("Number of hits in Photon Cluster").fill(pho.getCalorimeterHits().size());
+// aida.cloud1D("Energy of Photon Cluster").fill(pho.getEnergy());
+// aida.cloud2D("Num hits vs E of Photons").fill(pho.getEnergy(),pho.getCalorimeterHits().size());
+// }
+// if (phoclusters.size()==1) aida.cloud1D("Single Photon ESum").fill(PFAPhoE);
+// if (phoclusters.size()>0) aida.cloud1D("Total Cluster Photon ESum").fill(PFAPhoE);
+// }
+ } // end of event loop
+
+ private double[] layerEnergies(BasicCluster clus)
+ {
+ //FIXME could reuse this array
+ double[] layerEnergies = new double[_nLayers];
+ double clusterEnergy = 0.;
+ List<CalorimeterHit> hits = clus.getCalorimeterHits();
+// CalorimeterIDDecoder _decoderB = hit
+// CalorimeterIDDecoder _decoderB = (CalorimeterIDDecoder) event.getMetaData(chitsB).getIDDecoder();
+ //System.out.println("New cluster with "+hits.size()+ " hits"); // and energy "+clus.getEnergy());
+ for(CalorimeterHit hit : hits)
+ {
+ IDDecoder _decoder = hit.getIDDecoder();
+ _decoder.setID(hit.getCellID());
+ double e = hit.getRawEnergy();
+ int layer =_decoder.getLayer();
+// System.out.println("layer "+layer+" energy "+e);
+ clusterEnergy+=e;
+ layerEnergies[layer]+=e;
+ }
+// System.out.println("Layer 2 energy " + layerEnergies[2]);
+ //System.out.println("Cluster energy = " + clusterEnergy + "\n");
+
+ for(int i=0; i<layerEnergies.length; ++i)
+ {
+ layerEnergies[i]/=clusterEnergy;
+ }
+ //FIXME sum of cell energies does not equal cluster energy!
+// System.out.println(clusterEnergy+" "+clus.getEnergy());
+ return layerEnergies;
+ }
+
+ public void setInputClusterName(String clname)
+ {
+ _clusname = clname;
+ }
+ public void setOutputClusterName(String oclname)
+ {
+ _oclname = oclname;
+ }
+
+}
+
