lcsim/src/org/lcsim/contrib/SteveMagill

diff -N ChClusIDDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ChClusIDDriver.java	29 May 2008 21:33:13 -0000	1.1
@@ -0,0 +1,156 @@

+package org.lcsim.contrib.compile.SteveMagill;
+
+//  Driver to ID clusters based on MC input to the hits in the cluster - cluster ID is
+//  photon, neutral hadron, charged particle based on contributed energy
+
+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 ChClusIDDriver extends Driver
+{
+   private AIDA aida = AIDA.defaultInstance();
+   private ITree _tree;
+   private String _clusname;
+   private String _opclname;
+   private String _onclname;
+   private String _occlname;
+      
+   public ChClusIDDriver()
+   {
+
+   }
+   
+    protected void process(EventHeader event)
+    {
+        super.process(event);  // executes all added drivers
+            
+        List<BasicCluster> CHclusters = new ArrayList<BasicCluster>();  // CH clusters
+        List<BasicCluster> NHclusters = new ArrayList<BasicCluster>();  // NH clusters
+        List<BasicCluster> Phclusters = new ArrayList<BasicCluster>();  //  photons clusters
+
+        try
+        {
+            List<BasicCluster> inclusters = event.get(BasicCluster.class,_clusname);                
+
+            //  loop over all clusters
+            for (BasicCluster inclus : inclusters)
+            {
+                int numhits = inclus.getSize();
+                double clE = inclus.getEnergy();
+                double clphoE = 0.;
+                double clnhE = 0.;
+                double clchE = 0.;
+                double clphoh = 0.;
+                double clnhh = 0.;
+                double clchh = 0.;
+                List<CalorimeterHit> clhits = inclus.getCalorimeterHits();
+                for (CalorimeterHit clhit : clhits)
+                {
+                    SimCalorimeterHit mchit = (SimCalorimeterHit) clhit;
+                    double hitch = Math.abs(mchit.getMCParticle(0).getCharge());
+                    int hitid = mchit.getMCParticle(0).getPDGID();
+                    double hitE = mchit.getMCParticle(0).getEnergy();
+                    if (hitch==0)
+                    {
+                        if (hitid==22) 
+                        {
+                            clphoE += hitE;
+                            clphoh++;
+                        }
+                        if (hitid!=22) 
+                        {
+                            clnhE += hitE;
+                            clnhh++;
+                        }
+                    } else
+                    {
+                        clchE += hitE;
+                        clchh++;
+                    }
+                }
+                double phor = clphoh/numhits;
+                double nhr = clnhh/numhits;
+                double chr = clchh/numhits;
+                if (phor>nhr && phor>chr)
+                {
+                    aida.cloud1D("Photon Cluster Hit Ratio").fill(phor);
+                    Phclusters.add(inclus);
+                }
+                if (nhr>phor && nhr>chr)
+                {
+                    aida.cloud1D("Neutral Hadron Cluster Hit Ratio").fill(nhr);
+                    NHclusters.add(inclus);
+                }
+                if (chr>nhr && chr>phor)
+                {
+                    aida.cloud1D("Charged Hadron Cluster Hit Ratio").fill(chr);
+                    CHclusters.add(inclus);
+                }
+            
+            }  // end of cluster loop
+        }
+        catch (java.lang.IllegalArgumentException ex)
+        {
+            System.out.println("No clusters for ChCLusID finder");
+        }
+        if (Phclusters.size()>0) event.put(_opclname,Phclusters);
+        if (NHclusters.size()>0) event.put(_onclname,NHclusters);
+        if (CHclusters.size()>0) event.put(_occlname,CHclusters);
+        
+        // add both types of clusters here
+        List<BasicCluster> phoclusters = new ArrayList<BasicCluster>();
+        List<BasicCluster> hmphos = event.get(BasicCluster.class,"HMPhoClusters");
+        for (BasicCluster phcl : Phclusters)
+        {
+            phoclusters.add(phcl);
+        }
+        for (BasicCluster hmpho : hmphos)
+        {
+            phoclusters.add(hmpho);
+        }
+        event.put("PhotonClusters",phoclusters);
+            
+    }  // end of event loop
+    
+    
+    public void setInputClusterName(String clname)
+    {
+        _clusname = clname;
+    }
+    public void setOutputPhoClusterName(String opclname)
+    {
+        _opclname = opclname;
+    }
+    public void setOutputNHClusterName(String onclname)
+    {
+        _onclname = onclname;
+    }        
+        public void setOutputCHClusterName(String occlname)
+    {
+        _occlname = occlname;
+    }        
+}
+