lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon

diff -N FullE4JetFinder.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ FullE4JetFinder.java	28 Jun 2011 14:24:05 -0000	1.1
@@ -0,0 +1,156 @@

+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+import hep.physics.jet.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class FullE4JetFinder extends Driver
+{
+    String RPLname;
+    EventShape es;
+    String outN;
+    public FullE4JetFinder(String rpln)
+    {
+        RPLname = rpln;
+        es = new EventShape();
+        outN = RPLname+"FullE4Jet";
+    }
+    public FullE4JetFinder()
+    {
+        es = new EventShape();
+    }
+
+    public void setRPLname(String RPLname)
+    {
+        this.RPLname = RPLname;
+        outN = RPLname+"FullE4Jet";
+    }
+    protected void process(EventHeader event)
+    {
+        List<ReconstructedParticle> outlist = new ArrayList<ReconstructedParticle>();
+        List<ReconstructedParticle>[] dijet = new List[2];
+        List<ReconstructedParticle>[][] jet = new List[2][2];
+        dijet[0] = new ArrayList<ReconstructedParticle>();
+        dijet[1] = new ArrayList<ReconstructedParticle>();
+        double[] djE = new double[2];
+        Hep3Vector[] djP = new Hep3Vector[2];
+        djP[0] = new BasicHep3Vector();
+        djP[1] = new BasicHep3Vector();
+        jet[0][0] = new ArrayList<ReconstructedParticle>();
+        jet[0][1] = new ArrayList<ReconstructedParticle>();
+        jet[1][0] = new ArrayList<ReconstructedParticle>();
+        jet[1][1] = new ArrayList<ReconstructedParticle>();
+        double[][] jE = new double[2][2];
+        Hep3Vector[][] jP = new Hep3Vector[2][2];
+        jP[0][0] = new BasicHep3Vector();
+        jP[0][1] = new BasicHep3Vector();
+        jP[1][0] = new BasicHep3Vector();
+        jP[1][1] = new BasicHep3Vector();
+        List<ReconstructedParticle> rpl = event.get(ReconstructedParticle.class,RPLname);
+        Map<HepLorentzVector,ReconstructedParticle> ltormap = new HashMap<HepLorentzVector,ReconstructedParticle>();
+        for(ReconstructedParticle rp:rpl)
+        {
+            ltormap.put(new BasicHepLorentzVector(rp.getEnergy(),rp.getMomentum()),rp);
+        }
+        // Find thrust axis, and divide particles with hemisphere cut
+        es.setEvent(ltormap.keySet());
+        Hep3Vector thrusta = es.thrustAxis();
+        for(ReconstructedParticle rp:rpl)
+        {
+            int id = 1;
+            if(VecOp.dot(thrusta,rp.getMomentum()) > 0)id = 0;
+            dijet[id].add(rp);
+            djE[id] += rp.getEnergy();
+            djP[id] = VecOp.add(djP[id], rp.getMomentum());
+        }
+        // If less than 2 particles in a hemisphere, can't make 4 jets.
+        // Write empty list to event and return
+        if( (dijet[0].size() < 2)||(dijet[1].size() < 2) )
+        {
+            event.put(outN,outlist,ReconstructedParticle.class,0);
+            return;
+        }
+        // For each dijet, boost to cm and apply hemisphere cut there to
+        // divide into 2 jets.
+        for(int i=0;i<2;i++)
+        {
+            HepLorentzVector b4 = new BasicHepLorentzVector(djE[i],djP[i]);
+            Map<ReconstructedParticle,HepLorentzVector> rlmap = new HashMap<ReconstructedParticle,HepLorentzVector>();
+            for(ReconstructedParticle rp:dijet[i])
+            {
+                rlmap.put(rp,VecOp.boost(rp.asFourVector(),b4));
+            }
+            es.setEvent(rlmap.values());
+            Hep3Vector ta = es.thrustAxis();
+            for(ReconstructedParticle rp:dijet[i])
+            {
+                int id = 1;
+                if(VecOp.dot(rlmap.get(rp).v3(),ta) > 0)id = 0;
+                jet[i][id].add(rp);
+                jE[i][id] += rp.getEnergy();
+                jP[i][id] = VecOp.add(jP[i][id],rp.getMomentum());
+            }
+            if( (jet[i][0].size() < 2)||(jet[i][1].size() < 2) )
+            {
+                event.put(outN,outlist,ReconstructedParticle.class,0);
+                return;
+            }
+            if(jE[i][1] > jE[i][0])
+            {
+                double tE = jE[i][1];
+                Hep3Vector tP = jP[i][1];
+                List<ReconstructedParticle> tj = jet[i][0];
+                jE[i][1] = jE[i][0];
+                jE[i][0] = tE;
+                jP[i][1] = jP[i][0];
+                jP[i][0] = tP;
+                jet[i][1] = jet[i][0];
+                jet[i][0] = tj;
+            }
+        }
+        // Now have the 4 jets, make jet RP, sort, and write to event
+        BaseReconstructedParticle[] rjet = new BaseReconstructedParticle[4];
+        int ij = 0;
+        double maxE = 0.;
+        int maxij = -1;
+        for(int i=0;i<2;i++)
+        {
+            for(int j=0;j<2;j++)
+            {
+                rjet[ij] = new BaseReconstructedParticle(jE[i][j],jP[i][j]);
+                for(ReconstructedParticle rp:jet[i][j])
+                {
+                    rjet[ij].addParticle(rp);
+                }
+                if(jE[i][j] > maxE)
+                {
+                    maxE = jE[i][j];
+                    maxij = ij;
+                }
+                ij++;
+            }
+        }
+        int[] order = new int[4];
+        order[0] = maxij;
+        order[1] = maxij+1;
+        order[2] = (maxij+2)%4;
+        order[3] = order[2]+1;
+        for(int i=0;i<4;i++)
+        {
+            outlist.add(rjet[order[i]]);
+        }
+        event.put(outN,outlist,ReconstructedParticle.class,0);
+    }
+}