lcsim/src/org/lcsim/contrib/SteveMagill
diff -N ClMergeDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClMergeDriver.java 20 Aug 2008 20:23:55 -0000 1.1
@@ -0,0 +1,186 @@
+package org.lcsim.contrib.SteveMagill;
+
+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 ClMergeDriver extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ private double _clmeh;
+ private double _clmee;
+ private double _clmhh;
+ private int _minhits;
+ private ITree _tree;
+ private boolean mergD = false;
+ private String _emclusname;
+ private String _hadclusname;
+ private String _oclname;
+
+ public ClMergeDriver(double clmee, double clmhh, double clmeh, int minhits)
+ {
+ _clmeh = clmeh;
+ _clmee = clmee;
+ _clmhh = clmhh;
+ _minhits = minhits;
+ }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event); // executes all added drivers
+
+ // make a copies of em, had cluster lists since modifying and sort by E
+ List<BasicCluster> cEMclusters = new ArrayList<BasicCluster>();
+ List<BasicCluster> cHADclusters = new ArrayList<BasicCluster>();
+ try
+ {
+ List<BasicCluster> EMclusters = event.get(BasicCluster.class,_hadclusname);
+// System.out.println("Number of EM Clus" + EMclusters.size());
+ for (BasicCluster emclus : EMclusters)
+ {
+ cEMclusters.add(emclus);
+ }
+ Collections.sort(cEMclusters, new ClusterESort());
+ }
+ catch (java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No em clusters into ClMerge");
+ }
+
+ try
+ {
+ List<BasicCluster> HADclusters = event.get(BasicCluster.class,_hadclusname);
+// System.out.println("Number of HAD Clus" + HADclusters.size());
+ for (BasicCluster hadclus : HADclusters)
+ {
+ cHADclusters.add(hadclus);
+ }
+ Collections.sort(cHADclusters, new ClusterESort());
+ }
+ catch (java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No had clusters into ClMerge");
+ }
+
+ // loop over all EM clusters
+ List<BasicCluster> nhclusters = new ArrayList<BasicCluster>(); // neuh clusters
+ int nmerge = 0;
+ try
+ {
+ List<BasicCluster> EMclusters = event.get(BasicCluster.class,_emclusname);
+// System.out.println("Number of EM Clus" + EMclusters.size());
+ for (BasicCluster emclus : EMclusters)
+ {
+ BasicCluster newclus = new BasicCluster();
+ newclus.addCluster(emclus);
+ double ecp[] = emclus.getPosition();
+ double ecpx = ecp[0];
+ double ecpy = ecp[1];
+ double ecpz = ecp[2];
+ double ecR = Math.sqrt(ecpx*ecpx+ecpy*ecpy);
+ double eclth = Math.atan(ecR/ecpz);
+ if (eclth<0) eclth+=Math.PI;
+ double eclph = Math.atan2(ecpy,ecpx);
+ if (eclph<0) eclph+=2*Math.PI;
+ if (mergD) aida.cloud1D("Theta of NH Cand Clusters").fill(eclth);
+ if (mergD) aida.cloud1D("Phi of NH Cand Clusters").fill(eclph);
+
+ // now check with had clusters to make mergers
+ try
+ {
+// List<BasicCluster> HADclusters = event.get(BasicCluster.class,_hadclusname);
+// for (BasicCluster hadclus : HADclusters)
+ for (Iterator<BasicCluster> hcl = cHADclusters.iterator(); hcl.hasNext();)
+ {
+ BasicCluster ihcl = hcl.next();
+ double hcp[] = ihcl.getPosition();
+ double hcpx = hcp[0];
+ double hcpy = hcp[1];
+ double hcpz = hcp[2];
+ double hcR = Math.sqrt(hcpx*hcpx+hcpy*hcpy);
+ double hclth = Math.atan(hcR/hcpz);
+ if (hclth<0) hclth+=Math.PI;
+ double hclph = Math.atan2(hcpy,hcpx);
+ if (hclph<0) hclph+=2*Math.PI;
+ if (mergD) aida.cloud1D("Theta of NH Cand Clusters").fill(hclth);
+ if (mergD) aida.cloud1D("Phi of NH Cand Clusters").fill(hclph);
+ // calculate distance in theta phi
+ double dccth = Math.abs(eclth-hclth);
+ double dccph = Math.abs(eclph-hclph);
+ if (dccph>Math.PI) dccph = 2*Math.PI-dccph;
+ double dist = Math.sqrt(dccth*dccth+dccph*dccph);
+ if (mergD) aida.cloud1D("Distance between EM and HAD clusters").fill(dist);
+ // test for distance and merge if close, removing merged hadclusters from list
+ if (dist<_clmeh)
+ {
+ newclus.addCluster(ihcl);
+ hcl.remove();
+ nmerge++;
+ }
+ } // end of had cluster loop
+ }
+ catch (java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No had clusters in merger");
+ }
+
+ if (newclus.getSize()>_minhits) nhclusters.add(newclus);
+ } // end of EM cluster loop
+ }
+ catch (java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No em clusters in merger");
+ }
+
+ // now, add in remaining had clusters that were not merged
+ if (cHADclusters.size()>0)
+ {
+ for (BasicCluster rhcl : cHADclusters)
+ {
+ if (rhcl.getSize()>_minhits) nhclusters.add(rhcl);
+ }
+ }
+// System.out.println("Number of HAD merges" + nmerge);
+// System.out.println("Number of merged clusters" + nhclusters.size());
+
+ event.put(_oclname,nhclusters);
+
+ } // end of event loop
+
+ public void setInputEMClusterName(String emclname)
+ {
+ _emclusname = emclname;
+ }
+ public void setInputHADClusterName(String hadclname)
+ {
+ _hadclusname = hadclname;
+ }
+ public void setOutputClusterName(String oclname)
+ {
+ _oclname = oclname;
+ }
+
+}
+
