diff -N QqbarBarrelPhiAnalysisDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ QqbarBarrelPhiAnalysisDriver.java	5 Aug 2010 18:37:17 -0000	1.1
@@ -0,0 +1,221 @@

+package org.lcsim.contrib.Cassell.recon.analysis;
+import org.lcsim.recon.analysis.RMS90Calculator;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import hep.physics.vec.*;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.*;
+import java.util.*;
+import java.io.IOException;
+public class QqbarBarrelPhiAnalysisDriver extends Driver
+{
+   private AIDA aida = AIDA.defaultInstance();
+   int ievt;
+   int nmax = 1000000;
+   double[] phil = {3.333,6.667,10.,13.333,16.667,20.,23.333,26.667,30.};
+   double[] phiv;
+   double phidf;
+   double ctcut = 0.8;
+   double ZE;
+   RMS90Calculator calc;
+   IAnalysisFactory  analysisFactory;
+//   ITreeFactory      treeFactory;
+   ITree             tree;
+   IDataPointSetFactory dpsf ;
+   boolean first;
+   String pre;
+   String pre2;
+//
+// Keep track of the different cm Energies for plots
+//
+   int nE;
+   int[] Eevt;
+   public QqbarBarrelPhiAnalysisDriver()
+   {
+      ievt = 0;
+      analysisFactory = IAnalysisFactory.create();
+//      treeFactory = analysisFactory.createTreeFactory();
+//      tree = treeFactory.create();
+      tree = aida.tree();
+      dpsf   = analysisFactory.createDataPointSetFactory(tree);
+      calc = new RMS90Calculator();
+      first = true;
+      nE = 0;
+      Eevt = new int[10];
+      phiv = new double[phil.length];
+      double lv = 0.;
+      for(int i=0;i<phil.length;i++)
+      {
+         phiv[i] = lv + (phil[i] - lv)/2.;
+         lv = phil[i];
+      }
+   }
+   protected void process(EventHeader event)
+   {
+      List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+      double Zmass = 0.;
+      ZE = 0.;
+      double ct = 0.;
+      for(MCParticle p:mcl)
+      {
+         int id = Math.abs(p.getPDGID());
+         if( (id == 1 )||(id == 2 )||(id == 3 ) )
+         {
+            List<MCParticle> parents = p.getParents();
+            for(MCParticle parent:parents)
+            {
+               if(parent.getPDGID() == 23 )
+               {
+                   Zmass = parent.getMass();
+                   ZE = parent.getEnergy();
+                   Hep3Vector pp = p.getMomentum();
+                   ct = Math.abs(pp.z()/pp.magnitude());
+                   double phi = Math.atan2(pp.y(),pp.x());
+                   if(phi < 0.)phi += 2.*Math.PI;
+                   double phideg = phi*180./Math.PI;
+                   int n12 = (int) (phideg/30.);
+                   phidf = phideg - 30.*n12;
+               }
+            }
+         }
+      }
+//      System.out.println("Processing event "+ievt);
+      if(ct > ctcut)
+      {
+         ievt++;
+         return;
+      }
+      super.process(event);
+      int phib = 0;
+      for(int i=0;i<phil.length;i++)
+      {
+         if(phidf < phil[i])break;
+         phib++;
+      }
+      pre2 = "phibin "+phib+"/";
+//
+// Initialize cmE on first event
+//
+      if(first)
+      {
+         first = false;
+         int El = (int) (ZE + .5);
+         Eevt[0] = El;
+         nE = 1;
+//         tree.mkdirs("cmE="+Eevt[0]+"/");
+         pre = "cmE="+Eevt[0]+"/";
+         System.out.println("New CM Energy = "+El+"GeV at event "+ievt);
+      }
+//
+// See if cmE changed
+//
+      int Ec = (int) (ZE + .5);
+      boolean match = false;
+      for(int i=0;i<nE;i++)
+      {
+         if(Ec == Eevt[i])
+         {
+            match = true;
+            break;
+         }
+      }
+      if(!match)
+      {
+         System.out.println("New CM Energy = "+Ec+"GeV at event "+ievt);
+         Eevt[nE] = Ec;
+//         tree.mkdirs("cmE="+Eevt[nE]+"/");
+         pre = "cmE="+Eevt[nE]+"/";
+         nE++;
+      }
+      List<ReconstructedParticle> rl = event.get(ReconstructedParticle.class,"ReconstructedParticles");
+      double Etot = 0.;
+      double px = 0.;
+      double py = 0.;
+      double pz = 0.;
+      for(ReconstructedParticle p:rl)
+      {
+         Etot += p.getEnergy();
+         Hep3Vector mom = p.getMomentum();
+         px += mom.x();
+         py += mom.y();
+         pz += mom.z();
+      }
+      double M = Math.sqrt(Etot*Etot - px*px - py*py - pz*pz);
+      double dE = Etot - ZE;
+      double dM = M - Zmass;
+      aida.cloud1D(pre+"Event Energy",nmax).fill(Etot);
+      aida.cloud1D(pre+"Event Mass",nmax).fill(M);
+      aida.cloud1D(pre+"delta Event Energy",nmax).fill(dE);
+      aida.cloud1D(pre+"delta Event Mass",nmax).fill(dM);
+      aida.cloud1D(pre+pre2+"Event Energy",nmax).fill(Etot);
+      aida.cloud1D(pre+pre2+"Event Mass",nmax).fill(M);
+      aida.cloud1D(pre+pre2+"delta Event Energy",nmax).fill(dE);
+      aida.cloud1D(pre+pre2+"delta Event Mass",nmax).fill(dM);
+      ievt++;
+   }
+   protected void suspend()
+   {
+      IDataPointSet dpsE;
+      IDataPointSet[] dpsE2 = new IDataPointSet[nE];
+      IDataPointSet dpsM;
+      IDataPointSet[] dpsM2 = new IDataPointSet[nE];
+      dpsE = dpsf.create("dE over E vs cmE","dEoE vs E",2);
+      dpsM = dpsf.create("mean90 vs cmE","m90 vs E",2);
+      for(int iE=0;iE<nE;iE++)
+      {
+         dpsE2[iE] = dpsf.create("E="+Eevt[iE]+":dE over E vs phibin","Ecm="+Eevt[iE]+":dEoE vs phi",2);
+         dpsM2[iE] = dpsf.create("E="+Eevt[iE]+":mean90 vs phibin","Ecm="+Eevt[iE]+":m90 vs phi",2);
+         ZE = Eevt[iE];
+         String rname = "cmE="+Eevt[iE]+"/Event Energy";
+          ICloud1D thisc = aida.cloud1D(rname);
+          double r90 = calc.calculateRMS90(thisc);
+          double m90 = calc.getMEAN90();
+          int ent = thisc.entries();
+          double m = thisc.mean();
+          double r = thisc.rms();
+          double dEoE = r90/(m90);
+          double err2 = (1.1*r90/Math.sqrt(.9*ent));
+          double err = (1.1*r90/Math.sqrt(1.8*ent))/(m90);
+          dpsE.addPoint();
+          dpsM.addPoint();
+          IDataPoint dp = dpsE.point(iE);
+          IDataPoint dp2 = dpsM.point(iE);
+          dp.coordinate(0).setValue(ZE);
+          dp2.coordinate(0).setValue(ZE);
+          dp.coordinate(1).setValue(dEoE);
+          dp2.coordinate(1).setValue(m90);
+          dp.coordinate(1).setErrorPlus(err);
+          dp.coordinate(1).setErrorMinus(err);
+          dp2.coordinate(1).setErrorPlus(err2);
+          dp2.coordinate(1).setErrorMinus(err2);
+          for(int i=0;i<phil.length;i++)
+          {
+             rname = "cmE="+Eevt[iE]+"/phibin "+i+"/Event Energy";
+             thisc = aida.cloud1D(rname);
+             r90 = calc.calculateRMS90(thisc);
+             m90 = calc.getMEAN90();
+             ent = thisc.entries();
+             m = thisc.mean();
+             r = thisc.rms();
+             dEoE = r90/(m90);
+             err2 = (1.1*r90/Math.sqrt(.9*ent));
+             err = (1.1*r90/Math.sqrt(1.8*ent))/(m90);
+             dpsE2[iE].addPoint();
+             dpsM2[iE].addPoint();
+             dp = dpsE2[iE].point(i);
+             dp2 = dpsM2[iE].point(i);
+             dp.coordinate(0).setValue(phiv[i]);
+             dp2.coordinate(0).setValue(phiv[i]);
+             dp.coordinate(1).setValue(dEoE);
+             dp2.coordinate(1).setValue(m90);
+             dp.coordinate(1).setErrorPlus(err);
+             dp.coordinate(1).setErrorMinus(err);
+             dp2.coordinate(1).setErrorPlus(err2);
+             dp2.coordinate(1).setErrorMinus(err2);
+          }
+      }
+      super.suspend();
+   }
+}