lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon/analysis
diff -N QqbarBarrelPhiNonprojAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ QqbarBarrelPhiNonprojAnalysisDriver.java 10 Aug 2010 14:45:34 -0000 1.1
@@ -0,0 +1,222 @@
+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 QqbarBarrelPhiNonprojAnalysisDriver 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 QqbarBarrelPhiNonprojAnalysisDriver()
+ {
+ 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;
+ if(n12%2 == 1)phidf = 30 - phidf;
+ }
+ }
+ }
+ }
+// 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();
+ }
+}