diff -N PFASelect.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ PFASelect.java	18 May 2007 19:59:11 -0000	1.1
@@ -0,0 +1,170 @@

+package org.lcsim.contrib.SteveMagill;
+
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+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 PFASelect extends Driver
+{
+    
+    private AIDA aida = AIDA.defaultInstance();
+
+    public PFASelect()
+    {
+        add(new PFATemplate());
+        
+       //  now find jets with real Reco Particles
+       JetDriver PFAjet = new JetDriver();
+       PFAjet.setInputCollectionName("AllRecoParticles");
+       PFAjet.setOutputCollectionName("PFAJets");
+       JetFinder tworealjets = new FixNumberOfJetsFinder(2);
+       PFAjet.setFinder(tworealjets);
+       //  for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+//       JetFinder rktjet = new DurhamJetFinder(0.0005);
+//       PFAjet.setFinder(rktjet);
+       add(PFAjet);
+       
+       // make perfect jets out of perfect particles
+       JetDriver j = new JetDriver();
+       j.setInputCollectionName("PerfectRecoParticles");
+       j.setOutputCollectionName("PerfectJets");
+       JetFinder twojet = new FixNumberOfJetsFinder(2);
+       j.setFinder(twojet);
+       //  for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+//       JetFinder ktjet = new DurhamJetFinder(0.0005);
+//       j.setFinder(ktjet);
+       add(j);
+       
+       //  make plots etc. from results
+       PPFAPlotDriver ppfaplots = new PPFAPlotDriver();
+       ppfaplots.setMCPartName("CheatReconstructedParticles");
+       ppfaplots.setPerfJetName("PerfectJets");
+       ppfaplots.setPFAJetName("PFAJets");
+       ppfaplots.setPerfPartName("PerfectRecoParticles");
+       ppfaplots.setPFAPartName("AllRecoParticles");
+       add(ppfaplots);
+    }
+    
+    protected void process(EventHeader event)
+    {
+        //   Only process events with both jets costheta < costhetaMAX
+
+        List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+        boolean keepit = true;
+        double Zmass = 0.;
+        int ievt = 0;
+        double costmax = 0.8;
+        for(MCParticle p:mcl)
+        {
+           int id = Math.abs(p.getPDGID());
+           if( (id == 1 )||(id == 2 )||(id == 3 ) )
+           {
+              if(p.getParents().get(0).getPDGID() == 23 )
+              {
+                 Zmass = p.getParents().get(0).getMass();
+//                 System.out.println("Z Mass " + Zmass);
+                 Hep3Vector pp = p.getMomentum();
+                 double cost = Math.abs(pp.z()/pp.magnitude());
+                 if(cost > costmax)
+                 {
+                    keepit = false;
+                    continue;
+                 }
+              }
+           }
+        }
+        
+        if(!keepit)
+        {
+           ievt++;
+           return;
+        }
+        super.process(event);
+        //  add some plots here - get jet stuff
+        List<ReconstructedParticle> perfjets = event.get(ReconstructedParticle.class,"PerfectJets");
+        List<ReconstructedParticle> pfajets = event.get(ReconstructedParticle.class,"PFAJets");
+
+        //  here's some jet stuff for perfect jets
+        double jESum = 0;
+        double jPXSum = 0;
+        double jPYSum = 0;
+        double jPZSum = 0;
+        double djm = 0;
+        double[] jpz = new double[2];
+        double[] jp = new double[2];
+        int j = 0;
+        aida.cloud1D("Number of Perfect Jets").fill(perfjets.size());
+        if (perfjets.size() == 2)
+        {
+            //  calculate dijet mass - don't we have something that does this?
+            for (ReconstructedParticle perfjet : perfjets)
+            {   
+                aida.cloud1D("Perfect Jet Energy").fill(perfjet.getEnergy());
+                aida.cloud1D("Perfect Jet Mass").fill(perfjet.getMass());
+//                jESum += perfjet.getEnergy();
+                Hep3Vector j3v = perfjet.getMomentum();
+                jPXSum += j3v.x();
+                jPYSum += j3v.y();
+                jPZSum += j3v.z();
+                jp[j] = Math.sqrt(j3v.x()*j3v.x()+j3v.y()*j3v.y()+j3v.z()*j3v.z());
+                // force jet mass 0
+                jESum += jp[j];
+                aida.cloud1D("Perfect Jet Energy zero mass").fill(jp[j]);
+                jpz[j] = j3v.z();
+                 j++;
+            }    
+            double jE2 = jESum*jESum;
+            double jP2 = jPXSum*jPXSum+jPYSum*jPYSum+jPZSum*jPZSum;
+            djm += Math.sqrt(jE2-jP2);
+            aida.cloud1D("Perfect Dijet Mass").fill(djm);
+            aida.cloud1D("Diff Perfect Dijet Mass").fill((djm-Zmass)/Math.sqrt(Zmass)); 
+        }
+      
+        //  now real PFA jets
+        double rjESum = 0;
+        double rjPXSum = 0;
+        double rjPYSum = 0;
+        double rjPZSum = 0;
+        double rdjm = 0;
+        double[] jrpz = new double[2];
+        double[] jrp = new double[2];
+        int k = 0;
+        aida.cloud1D("Number of PFA Jets").fill(pfajets.size());
+        if (pfajets.size() == 2)
+        {
+            //  calculate dijet mass - don't we have something that does this?
+            for (ReconstructedParticle pfajet : pfajets)
+            {
+                aida.cloud1D("PFA Jet Energy").fill(pfajet.getEnergy());
+                aida.cloud1D("PFA Jet Mass").fill(pfajet.getMass());
+//                rjESum += pfajet.getEnergy();
+                Hep3Vector rj3v = pfajet.getMomentum();
+                rjPXSum += rj3v.x();
+                rjPYSum += rj3v.y();
+                rjPZSum += rj3v.z();
+                jrp[k] = Math.sqrt(rj3v.x()*rj3v.x()+rj3v.y()*rj3v.y()+rj3v.z()*rj3v.z());
+                // force jet mass 0
+                rjESum += jrp[k];
+                aida.cloud1D("PFA Jet Energy zero mass").fill(jrp[k]);
+                jrpz[k] = rj3v.z();
+                k++;
+            }
+            double rjE2 = rjESum*rjESum;
+            double rjP2 = rjPXSum*rjPXSum+rjPYSum*rjPYSum+rjPZSum*rjPZSum;
+            rdjm += Math.sqrt(rjE2-rjP2);
+            aida.cloud1D("PFA Dijet Mass").fill(rdjm);
+            aida.cloud1D("Diff PFA Dijet Mass").fill((rdjm-Zmass)/Math.sqrt(Zmass));
+            aida.cloud1D("Difference PFAJet and PerfJet DiJet Mass").fill(rdjm-djm);
+        }
+    }
+}