diff -N MCHiggs.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MCHiggs.java	19 Dec 2012 17:14:06 -0000	1.1
@@ -0,0 +1,81 @@

+package org.lcsim.mcd.analysis;
+
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+/**
+ *
+ * @author agias
+ */
+public class MCHiggs extends Driver{
+    
+    public MCHiggs()
+    {
+    }
+    
+    private AIDA aida = AIDA.defaultInstance();
+    
+     protected String mcParticlesSkimmedName = "MCParticle";
+
+    public void process(EventHeader event) {
+	List<MCParticle> particles = event.get(MCParticle.class,mcParticlesSkimmedName);
+        List<MCParticle> children = new ArrayList<MCParticle>();
+        
+       boolean hfound = false;
+	for (MCParticle particle : particles) {
+            // 25 == h0/H01
+	    if ((particle.getPDGID()==25)&&(!hfound)) {
+                children.addAll(particle.getDaughters());
+                hfound=true;
+            }
+	}
+        
+        for (MCParticle particle : children) {
+            /*
+             * Only count positive PDGID's because negatives 
+             * are just antiparticles
+             * 4    == c
+             * 5    == b
+             * 15   == tau-
+             * 21   == gluon
+             * 23   == Zo
+             * 24   == W
+             */
+            if (particle.getPDGID() > 0) {
+                aida.cloud1D("decay modes").fill(particle.getPDGID());
+            }
+            if (particle.getPDGID() == 5) {
+                Hep3Vector p = particle.getMomentum();
+                Hep3Vector p_bar = VecOp.neg(p);
+                
+                double mass = particle.getMass();
+                double en   = particle.getEnergy();
+                
+                Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
+                                
+                double p_m = p.magnitude();
+                        
+                double pt = VecOp.dot(p, z_unit);
+                double pt_bar = VecOp.dot(p_bar, z_unit);
+                
+                aida.cloud1D("b momentum").fill(p_m);
+                aida.cloud1D("b transverse momentum").fill(pt);
+                aida.cloud1D("b_bar transverse momentum").fill(pt_bar);
+                
+                // b and b_bar are produced with equal energy, opposite momenta
+                double inv_mass = Math.sqrt(2.*mass*mass + 2.*(en*en-VecOp.dot(p, p_bar)));
+                aida.cloud1D("H invariant mass").fill(inv_mass);
+            }
+        }
+        aida.cloud1D("num children").fill(children.size());
+    }
+    
+}