diff -N HiggsRecoilMassAnalysis.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ HiggsRecoilMassAnalysis.java	30 Apr 2013 23:07:19 -0000	1.1
@@ -0,0 +1,101 @@

+package org.lcsim.analysis;
+
+import hep.aida.ITree;
+import hep.physics.vec.BasicHepLorentzVector;
+import hep.physics.vec.HepLorentzVector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+import static java.lang.Math.abs;
+
+/**
+ *
+ * @author Norman A Graf
+ * 
+ *  @version $Id:
+ */
+public class HiggsRecoilMassAnalysis extends Driver
+{
+    private AIDA aida = AIDA.defaultInstance();
+    private ITree _tree;
+    
+    // analysis-specific variables
+    private String _reconstructedParticleListName = "PandoraPFOCollection";
+    private double _minMuonMomentum = 20.;
+    private double _cmsEnergy = 250.;
+    /** Creates a new instance of HiggsRecoilAnalysis */
+    public HiggsRecoilMassAnalysis()
+    {
+    }
+    
+    protected void process(EventHeader event)
+    {
+//       Fetch the list of ReconstructedParticles...
+         List<ReconstructedParticle> rpList = event.get(ReconstructedParticle.class, _reconstructedParticleListName);      
+//       Create a collection to hold muons found in this event
+        List<ReconstructedParticle> muons = new ArrayList<ReconstructedParticle>();
+        
+        for (ReconstructedParticle rp : rpList)
+        {
+            if(abs(rp.getType())==13)
+            {
+                if(rp.getMomentum().magnitude()>_minMuonMomentum)
+                {
+                    muons.add(rp);
+                    aida.cloud1D("Muon Momentum").fill(rp.getMomentum().magnitude());
+                }
+            }
+        }
+//     Calculate the invariant mass of the two muons
+       double invMass = 0.;
+//        
+//     create the Z
+       HepLorentzVector Z = new  BasicHepLorentzVector();
+//        
+        aida.cloud1D("number of muons passing cut").fill(muons.size());
+//        
+        if(muons.size()==2)
+        {
+            for(int i=0; i<2; ++i)
+            {
+                ReconstructedParticle mu = muons.get(i);
+                Z = VecOp.add(Z, mu.asFourVector());
+            }
+            
+            invMass = Z.magnitude();
+            aida.cloud1D("Dimuon invariant mass").fill(invMass);
+            
+            // calculate the recoil mass if we have a Z
+            double zMass = 91.;
+            double zMassWindow = 5.;
+            if(Math.abs(zMass-invMass)<zMassWindow)
+            {
+                // M_recoil = sqrt( s^2 - 2*sqrt(s)*E_ll + M^2_ll )                
+                double s = _cmsEnergy;
+                double recoilMass = Math.sqrt(s*s -(2.*s*Z.t()) + Z.magnitudeSquared());
+                aida.cloud1D("Recoil mass").fill(recoilMass);
+                aida.histogram1D("Recoil mass binned and weighted", 100, 100.,150.).fill(recoilMass,event.getWeight());
+            }
+        }
+    }
+    
+    public void setReconstructedParticleListName(String rplname)
+    {
+        _reconstructedParticleListName = rplname;
+    }
+    
+    public void setMinMuonMomentum( double minMomentum)
+    {
+        _minMuonMomentum = minMomentum;
+    }
+    public void setCMSEnergy(double cmsEnergy)
+    {
+        _cmsEnergy = cmsEnergy;
+    }
+}
+