diff -N MCqqbarME.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MCqqbarME.java	12 Aug 2009 18:29:28 -0000	1.1
@@ -0,0 +1,125 @@

+package org.lcsim.contrib.SteveMagill;
+
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.util.Driver;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.jet.*;
+
+/*
+ * PFASelect includes fiducial cut to select region for PFA analysis
+ */
+
+public class MCqqbarME extends Driver
+{
+    private String _mcname;
+    private String _rpname;
+    private AIDA aida = AIDA.defaultInstance();
+
+    public MCqqbarME()
+    {
+        
+    }
+    
+    protected void process(EventHeader event)
+    {
+        super.process(event);
+
+        //  get true qqbar E and M and get Nayeli's mass correction
+        List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+        double qqbarM = 0.;
+        double qqbarE = 0.;
+        double mcpx = 0.;
+        double mcpy = 0.;
+        double mcpz = 0.;
+        double Etot = 0.;
+        double pxprime = 0.;
+        double pyprime = 0.;
+        double pzprime = 0.;
+        for(MCParticle p:mcl)
+        {
+           int gs = p.getGeneratorStatus();
+           int pid = Math.abs(p.getPDGID());
+           double pch = p.getCharge();
+//           if (gs==1 || (gs==0 && pid==2112))  //  only want neutrons which are called FS, so don't use this
+           if (gs==1 && (pid != 12 && pid != 14 && pid != 16))  //  all FS particles except neutrinos
+           {
+               Etot += p.getEnergy();
+               mcpx += p.getPX();
+               mcpy += p.getPY();
+               mcpz += p.getPZ();
+               Hep3Vector rv = p.getEndPoint();
+               double rmag = rv.magnitude();
+               pxprime += p.getEnergy()*rv.x()/rmag;
+               pyprime += p.getEnergy()*rv.y()/rmag;
+               pzprime += p.getEnergy()*rv.z()/rmag;
+           }
+           int id = Math.abs(p.getPDGID());
+           if( (id == 1 )||(id == 2 )||(id == 3 ) )
+           {
+              int parid = Math.abs(p.getParents().get(0).getPDGID());
+              if(parid == 24 || parid == 23)
+              {
+                 qqbarM = p.getParents().get(0).getMass();
+                 qqbarE = p.getParents().get(0).getEnergy();
+//                 System.out.println("qqbar Mass " + qqbarM);
+              }
+           }
+        }
+        double pptot = Math.sqrt(pxprime*pxprime+pyprime*pyprime+pzprime*pzprime);
+        double mprime = Math.sqrt(Etot*Etot-pptot*pptot);
+        double ptot = Math.sqrt(mcpx*mcpx+mcpy*mcpy+mcpz*mcpz);
+        double minv = Math.sqrt(Etot*Etot-ptot*ptot);
+        double delm = mprime-minv;
+        aida.cloud1D("MC FS Invariant Mass").fill(minv);
+        aida.cloud1D("MC FS Prime Invariant Mass").fill(mprime);
+        aida.cloud1D("Delta M Mp M").fill(delm);
+//        System.out.println("MC FS Invariant Mass " + minv);
+
+        aida.cloud1D("qqbar Parent Mass").fill(qqbarM);
+        aida.cloud1D("qqbar Parent Energy").fill(qqbarE);
+
+//  make reconstructed particles out of MC FS list
+        double epho = 0.;
+        double ETot = 0.;
+        List<ReconstructedParticle> rpList = new ArrayList<ReconstructedParticle>();
+        List<MCParticle> mcfsp = event.get(MCParticle.class,_mcname);
+        for (MCParticle mcf : mcfsp)
+        {
+            int mcid = Math.abs(mcf.getPDGID());
+//            System.out.println("PDGID " + mcid);
+            if (mcid != 12 && mcid != 14 && mcid != 16)
+            {
+                double mcE = mcf.getEnergy();
+                ETot += mcE;
+                if (mcid == 22) epho += mcE;
+                Hep3Vector mcP = mcf.getMomentum();
+                BaseReconstructedParticle rp = new BaseReconstructedParticle(mcE, mcP);
+                rp.setMass(mcf.getMass());
+                rp.setCharge(mcf.getCharge());
+                rpList.add(rp);
+            }
+        }
+        event.put(_rpname,rpList, ReconstructedParticle.class,0);
+        double phorat = epho/ETot;
+        aida.cloud1D("Fraction of photons").fill(phorat);
+
+    }
+
+  public void setMCParticleName(String mcname)
+  {
+      _mcname = mcname;
+  }
+
+  public void setMCRecoParticleName(String rpname)
+  {
+      _rpname = rpname;
+  }
+
+}