diff -u -r1.9 -r1.10
--- MCHiggs.java	19 Dec 2012 18:29:04 -0000	1.9
+++ MCHiggs.java	8 Feb 2013 16:39:52 -0000	1.10
@@ -10,6 +10,7 @@

 import hep.physics.vec.Hep3Vector;
 import hep.physics.vec.BasicHep3Vector;
 import hep.physics.vec.VecOp;

+import org.lcsim.mc.fast.MCFast;

 
 /**
  *

@@ -19,6 +20,7 @@

     
     public MCHiggs()
     {

+        add(new MCFast());

     }
     
     private AIDA aida = AIDA.defaultInstance();

@@ -37,6 +39,8 @@

                 hfound=true;
             }
 	}

+        MCParticle p1 = children.get(0);
+        MCParticle p2 = children.get(0);

         
         for (MCParticle particle : children) {
             /*

@@ -53,58 +57,76 @@

              */
             if ((particle.getPDGID() > 0) && (particle.getPDGID() != 25)) {
                 aida.cloud1D("decay modes").fill(particle.getPDGID());

+                p1 = particle;

             }

-            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.histogram1D("b momentum hist",100,62.3108,62.3126).fill(p_m);
-                //aida.cloud1D("b transverse momentum").fill(pt);
-                aida.histogram1D("b transverse momentum hist", 100,0,63).fill(pt);
-                //aida.cloud1D("b_bar transverse momentum").fill(pt_bar);
-                aida.histogram1D("b_bar transverse momentum hist", 100,0,63).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.histogram1D("H inv mass hist",100,124.990,124.995).fill(inv_mass);
-                aida.histogram1D("H inv mass hist big",200,100,130).fill(inv_mass);
-                
-                aida.cloud1D("b phi").fill(VecOp.phi(p));
-                aida.cloud1D("b cos theta").fill(VecOp.cosTheta(p));
-            }
-            if (particle.getPDGID() == 22) {
-                double mass = particle.getMass();
-                double ener = particle.getEnergy();
-                Hep3Vector p = particle.getMomentum();
-                Hep3Vector p2 = VecOp.neg(p);
-                double p_m = p.magnitude();
-                double p2_m = p_m;
-                
-                Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
-                
-                double pt = VecOp.dot(p, z_unit);
-                double pt2 = VecOp.dot(p2,z_unit);
-                
-                
-                double inv_mass = Math.sqrt(2*(ener*ener - VecOp.dot(p, p2)));
-                
-                aida.cloud1D("gamma momentum").fill(p_m);
-                aida.cloud1D("gamma transverse momentum").fill(pt);
-                aida.cloud1D("gamma higgs inv mass").fill(inv_mass);

+            if ((particle.getPDGID() < 0)) {
+                p2 = particle;

             }
         }

+        
+        if (p1.getPDGID() == 5) {
+            Hep3Vector p = p1.getMomentum();
+            Hep3Vector p_bar = p2.getMomentum();
+
+            double mass1 = p1.getMass();
+            double en1   = p1.getEnergy();
+            double mass2 = p2.getMass();
+            double en2   = p2.getEnergy();
+
+            Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
+
+            double p1_m = p.magnitude();
+            double p2_m = p_bar.magnitude();
+
+            double pt = VecOp.dot(p, z_unit);
+            double pt_bar = VecOp.dot(p_bar, z_unit);
+
+            aida.cloud1D("b momentum").fill(p1_m);
+            aida.cloud1D("b_bar momentum").fill(p2_m);
+            aida.histogram1D("b momentum hist",100,62.3108,62.3126).fill(p1_m);
+            aida.histogram1D("b_bar momentum hist",100,62.3,62.4).fill(p2_m);
+            
+            aida.cloud1D("b transverse momentum").fill(pt);
+            aida.cloud1D("b_bar transverse momentum").fill(pt_bar);
+            aida.histogram1D("b transverse momentum hist", 100,0,63).fill(pt);
+            aida.histogram1D("b_bar transverse momentum hist", 100,0,63).fill(pt_bar);
+
+
+            double inv_mass = Math.sqrt(mass1*mass1 + mass2*mass2 + 2.*(en1*en2-VecOp.dot(p, p_bar)));
+            aida.cloud1D("H invariant mass").fill(inv_mass);
+            aida.histogram1D("H inv mass hist",100,124.990,124.995).fill(inv_mass);
+            aida.histogram1D("H inv mass hist big",200,100,130).fill(inv_mass);
+
+            aida.cloud1D("b phi").fill(VecOp.phi(p));
+            aida.cloud1D("bbar phi").fill(VecOp.phi(p_bar));
+            aida.cloud1D("b cos theta").fill(VecOp.cosTheta(p));
+            aida.cloud1D("bbar cos theta").fill(VecOp.cosTheta(p_bar));
+        }
+        if (p1.getPDGID() == 22) {
+            double mass1 = p1.getMass();
+            double ener1 = p1.getEnergy();
+            double mass2 = p2.getMass();
+            double ener2 = p2.getEnergy();
+
+            Hep3Vector pt1 = p1.getMomentum();
+            Hep3Vector pt2 = p2.getMomentum();
+            double p1_m = pt1.magnitude();
+            double p2_m = pt2.magnitude();
+
+            Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
+
+            double p1trans = VecOp.dot(pt1, z_unit);
+            double p2trans = VecOp.dot(pt2,z_unit);
+
+
+            double inv_mass = Math.sqrt(2*(ener1*ener2 - VecOp.dot(pt1, pt2)));
+
+            aida.cloud1D("gamma1 momentum").fill(p1_m);
+            aida.cloud1D("gamma2 momentum").fill(p2_m);
+            aida.cloud1D("gamma1 transverse momentum").fill(p1trans);
+            aida.cloud1D("gamma2 transverse momentum").fill(p2trans);
+            aida.cloud1D("gamma higgs inv mass").fill(inv_mass);
+        }

         aida.cloud1D("num children").fill(children.size());
     }