diff -N MakePerfectPDG94JetsCT.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MakePerfectPDG94JetsCT.java	18 Jul 2011 17:11:20 -0000	1.1
@@ -0,0 +1,181 @@

+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.event.util.CreateFinalStateMCParticleList;
+
+/**
+ *
+ * @author cassell
+ */
+public class MakePerfectPDG94JetsCT extends Driver
+{
+    List<ReconstructedParticle> pjlo;
+    RPESort rpes;
+    public MakePerfectPDG94JetsCT()
+    {
+        rpes = new RPESort();
+        add(new CreateFinalStateMCParticleList("Gen"));
+    }
+    protected void process(EventHeader event)
+    {
+        super.process(event);
+        List<MCParticle> all = event.get(MCParticle.class,"MCParticle");
+        List<MCParticle> fs = event.get(MCParticle.class,"GenFinalStateParticles");
+        // Ordered list of perfect jets
+        pjlo = new ArrayList<ReconstructedParticle>();
+        // Create the perfect pdg94 ReconstructedParticles
+        List<ReconstructedParticle> p94l = new ArrayList<ReconstructedParticle>();
+        for(MCParticle p:all)
+        {
+            if(p.getPDGID() == 94)
+            {
+                // Make a ReconstructedParticle from the pdg94
+                ReconstructedParticle r94 = makeRP94(p,fs);
+                // Add it to the list
+                p94l.add(r94);
+                // Sort the jets for this pdg94 by energy
+                Collections.sort(r94.getParticles(),rpes);
+            }
+        }
+        event.put("PerfectPDG94", p94l);
+        // Order the jets: jet 0 = maxE jet, jet1 = jet0 partner,
+        // jet2 = maxE jet from other pdg94, jet3 = jet2 partner
+        List<ReconstructedParticle> t94l = new ArrayList<ReconstructedParticle>(p94l);
+        while(t94l.size() > 0)
+        {
+            double maxct = -1.;
+            int n = t94l.size();
+            int i94 = -1;
+            for(int i=0;i<n;i++)
+            {
+                double ct = t94l.get(i).getMomentum().z()/t94l.get(i).getMomentum().magnitude();
+                if(ct > maxct)
+                {
+                    maxct = ct;
+                    i94 = i;
+                }
+            }
+            for(ReconstructedParticle jet:t94l.get(i94).getParticles())
+            {
+                pjlo.add(jet);
+            }
+            t94l.remove(i94);
+        }
+        event.put("PerfectJets",pjlo,ReconstructedParticle.class,0);
+        List<ReconstructedParticle> pjpl = new ArrayList<ReconstructedParticle>();
+        for(ReconstructedParticle j:pjlo)
+        {
+            for(ReconstructedParticle rp:j.getParticles())
+            {
+                pjpl.add(rp);
+            }
+        }
+        event.put("PerfectJetParticles",pjpl,ReconstructedParticle.class,0);
+    }
+    public ReconstructedParticle makeRP94(MCParticle pdg94,List<MCParticle> fsl)
+    {
+        // Make the original ReconstructedParticle
+        BaseReconstructedParticle pdg94rp = rpFromMC(pdg94);
+        // Construct and add the quark ReconstructedParticles
+        int type = 0;
+        int nused = 0;
+        for(MCParticle p:pdg94.getParents())
+        {
+            if(pdg94.getParents().size() > 2)
+            {
+                if(p.getDaughters().size() > 1)continue;
+            }
+            // Make the original quark reconstructedParticle
+            BaseReconstructedParticle rp = rpFromMC(p);
+            // Add it to the pdg94 particle
+            pdg94rp.addParticle(rp);
+            rp.setType(p.getPDGID());
+            type += p.getPDGID();
+            nused++;
+        }
+        // Adjust the pdg94 type
+        if(nused != 2)type -= 10;
+        else
+        {
+            if(pdg94.getParents().size() > 2)
+            {
+                type += 10;
+            }
+        }
+        pdg94rp.setType(type);
+        // Now add the fs particles traced to this pdg94 to the appropriate quark
+        for(MCParticle p:fsl)
+        {
+            if(trace(p,pdg94))
+            {
+                ReconstructedParticle rp = rpFromMC(p);
+                findNearestQ(rp,pdg94rp.getParticles(),pdg94rp).getParticles().add(rp);
+            }
+        }
+        // Now replace the quark 4-Vector with the actual sum of its components
+        for(ReconstructedParticle q:pdg94rp.getParticles())
+        {
+            HepLorentzVector v4 = new BasicHepLorentzVector();
+            for(ReconstructedParticle p:q.getParticles())
+            {
+                v4 = VecOp.add(v4,p.asFourVector());
+            }
+            BaseReconstructedParticle bq = (BaseReconstructedParticle) q;
+            bq.set4Vector(v4);
+            bq.setMass(v4.magnitude());
+        }
+        return pdg94rp;
+    }
+    public BaseReconstructedParticle rpFromMC(MCParticle p)
+    {
+        BaseReconstructedParticle rp = new BaseReconstructedParticle(p.getEnergy(),p.getMomentum());
+        rp.setMass(p.getMass());
+        rp.setType(p.getPDGID());
+        double charge = 0.;
+        try
+        {
+            charge = p.getCharge();
+        }
+        catch (hep.physics.particle.properties.UnknownParticleIDException e)
+        {
+        }
+        rp.setCharge(charge);
+        return rp;
+    }
+    public boolean trace(MCParticle fsp, MCParticle p94)
+    {
+        MCParticle pp = fsp;
+        while(pp.getParents().size() > 0)
+        {
+            if(pp.getParents().get(0) == p94)return true;
+            pp = pp.getParents().get(0);
+        }
+        return false;
+    }
+    public ReconstructedParticle findNearestQ(ReconstructedParticle p, List<ReconstructedParticle> ql, ReconstructedParticle bp)
+    {
+        ReconstructedParticle np = null;
+        double maxct = -2.;
+        HepLorentzVector pl = VecOp.boost(p.asFourVector(),bp.asFourVector());
+        for(ReconstructedParticle q:ql)
+        {
+            double ct = VecOp.dot(VecOp.boost(q.asFourVector(),bp.asFourVector()).v3(),pl.v3());
+            if(ct > maxct)
+            {
+                maxct = ct;
+                np = q;
+            }
+        }
+        return np;
+    }
+}