lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon/analysis
diff -N QqbarAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ QqbarAnalysisDriver.java 4 Nov 2010 15:12:56 -0000 1.1
@@ -0,0 +1,351 @@
+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.*;
+public class QqbarAnalysisDriver extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ int ievt;
+ int nmax = 1000000;
+ String skip = "ForVtxingReconParticles";
+ List<String> rplnames;
+ double[] ctl = {.1,.2,.3,.4,.5,.6,.7,.8,.9,.925,.95,.975,1.};
+ String[] ctname;
+ double[] cscl = {1.005,1.02,1.05,1.1,1.2,1.4,1.67};
+ String[] cscname;
+ double[] cutval = {.8,.95,.95};
+ String[] cutname;
+ double ZE;
+ RMS90Calculator calc;
+ IAnalysisFactory analysisFactory;
+ ITreeFactory treeFactory;
+ ITree tree;
+ IDataPointSetFactory dpsf ;
+ boolean first;
+//
+// Keep track of the different cm Energies for plots
+//
+ int nE;
+ int[] Eevt;
+ public QqbarAnalysisDriver()
+ {
+ ievt = 0;
+ analysisFactory = IAnalysisFactory.create();
+ treeFactory = analysisFactory.createTreeFactory();
+// tree = treeFactory.create("DPS.aida","xml",false,true);
+ tree = treeFactory.create();
+ dpsf = analysisFactory.createDataPointSetFactory(tree);
+ calc = new RMS90Calculator();
+ rplnames = new ArrayList<String>();
+ ctname = new String[ctl.length];
+ cscname = new String[cscl.length];
+ cutname = new String[cutval.length];
+ double ll = 0.;
+ for(int i=0;i<ctl.length;i++)
+ {
+ if(i > 0)ll = ctl[i-1];
+ double ul = ctl[i];
+ double lls = ((int)(ll*100. + .5))/100.;
+ double uls = ((int)(ul*100. + .5))/100.;
+ ctname[i] = lls+"<ct<"+uls;
+ }
+ ll = 0.;
+ for(int i=0;i<cscl.length;i++)
+ {
+ if(i > 0)ll = cscl[i-1];
+ double ul = cscl[i];
+ double lls = ((int)(ll*1000. + .5))/1000.;
+ double uls = ((int)(ul*1000. + .5))/1000.;
+ cscname[i] = lls+"<csc<"+uls;
+ }
+ cutname[0] = "barrel";
+ cutname[1] = "forward";
+ cutname[2] = "all";
+ first = true;
+ nE = 0;
+ Eevt = new int[10];
+ }
+ 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 ) )
+ {
+ if(p.getParents().get(0).getPDGID() == 23 )
+ {
+ Zmass = p.getParents().get(0).getMass();
+ ZE = p.getParents().get(0).getEnergy();
+ Hep3Vector pp = p.getMomentum();
+ ct = Math.abs(pp.z()/pp.magnitude());
+ }
+ }
+ }
+// System.out.println("Processing event "+ievt);
+ super.process(event);
+//
+// Initialize cmE on first event
+//
+ if(first)
+ {
+ first = false;
+ int El = (int) (ZE + .5);
+ Eevt[0] = El;
+ nE = 1;
+ tree.mkdirs("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]+"/");
+ nE++;
+ }
+ int ctb = 0;
+ for(int i=0;i<ctl.length;i++)
+ {
+ if(ct < ctl[i])break;
+ ctb++;
+ }
+ double csct = 1./Math.sqrt(1. - ct*ct);
+ int cscb = 0;
+ for(int i=0;i<cscl.length;i++)
+ {
+ if(csct < cscl[i])break;
+ cscb++;
+ }
+ List<List<ReconstructedParticle>> all = event.get(ReconstructedParticle.class);
+ for(List<ReconstructedParticle> rl:all)
+ {
+ if(event.getMetaData(rl).getName().compareTo(skip) != 0)
+ {
+ String rpname = event.getMetaData(rl).getName();
+ if(!rplnames.contains(rpname))rplnames.add(rpname);
+ String rname = "cmE="+Ec+"/"+rpname;
+ double Etot = 0.;
+ double px = 0.;
+ double py = 0.;
+ double pz = 0.;
+ double Ech = 0.;
+ double Enh = 0.;
+ double Eph = 0.;
+ for(ReconstructedParticle p:rl)
+ {
+ Etot += p.getEnergy();
+ Hep3Vector mom = p.getMomentum();
+ px += mom.x();
+ py += mom.y();
+ pz += mom.z();
+ if(p.getCharge() == 0.)
+ {
+ if(p.getMass() == 0.)Eph += p.getEnergy();
+ else Enh += p.getEnergy();
+ }
+ else Ech += p.getEnergy();
+ }
+ double M = Math.sqrt(Etot*Etot - px*px - py*py - pz*pz);
+ double dE = Etot - ZE;
+/*
+ if(Etot/ZE > 1.5)
+ {
+ System.out.println("dE = "+dE+": Etot, ZE = "+Etot+", "+ZE);
+ throw new RuntimeException("oops");
+ }
+*/
+ double nhfrac = Enh/ZE;
+ double dM = M - Zmass;
+ int nhfb = 0;
+ if(nhfrac > .25)nhfb = 4;
+ else if(nhfrac > .1)nhfb = 3;
+ else if(nhfrac > .05)nhfb = 2;
+ else if(nhfrac > 0.)nhfb = 1;
+ if(ct < cutval[0])
+ {
+ aida.cloud1D(rname+"/"+cutname[0]+"/dE",nmax).fill(dE);
+ aida.cloud1D(rname+"/"+cutname[0]+"/dM",nmax).fill(dM);
+ aida.cloud1D(rname+"/"+cutname[0]+"/nhE fraction",nmax).fill(nhfrac);
+ aida.cloud1D(rname+"/"+cutname[0]+"/nhfbin "+nhfb+"/dE",nmax).fill(dE);
+ }
+ else if(ct < cutval[1])
+ {
+ aida.cloud1D(rname+"/"+cutname[1]+"/dE",nmax).fill(dE);
+ aida.cloud1D(rname+"/"+cutname[1]+"/dM",nmax).fill(dM);
+ aida.cloud1D(rname+"/"+cutname[1]+"/nhE fraction",nmax).fill(nhfrac);
+ aida.cloud1D(rname+"/"+cutname[1]+"/nhfbin "+nhfb+"/dE",nmax).fill(dE);
+ }
+ if(ct < cutval[2])
+ {
+ aida.cloud1D(rname+"/"+cutname[2]+"/dE",nmax).fill(dE);
+ aida.cloud1D(rname+"/"+cutname[2]+"/dM",nmax).fill(dM);
+ aida.cloud1D(rname+"/"+cutname[2]+"/nhE fraction",nmax).fill(nhfrac);
+ aida.cloud1D(rname+"/"+cutname[2]+"/nhfbin "+nhfb+"/dE",nmax).fill(dE);
+ }
+ aida.cloud1D(rname+"/"+ctname[ctb]+"/dE",nmax).fill(dE);
+ aida.cloud1D(rname+"/"+ctname[ctb]+"/dM",nmax).fill(dM);
+ if(ct < .8)
+ {
+ aida.cloud1D(rname+"/"+cscname[cscb]+"/dE",nmax).fill(dE);
+ aida.cloud1D(rname+"/"+cscname[cscb]+"/dM",nmax).fill(dM);
+ }
+ }
+ }
+ ievt++;
+ }
+ protected void suspend()
+ {
+ IDataPointSet[][] dpsE = new IDataPointSet[rplnames.size()][cutname.length];
+ IDataPointSet[][] dpsE2 = new IDataPointSet[rplnames.size()][cutname.length];
+ IDataPointSet[][] dpsE3 = new IDataPointSet[rplnames.size()][cutname.length];
+ for(int i=0;i<rplnames.size();i++)
+ {
+ for(int j=0;j<cutname.length;j++)
+ {
+ dpsE[i][j] = dpsf.create(rplnames.get(i)+":"+cutname[j]+":alpha90 vs jetE",rplnames.get(i)+":"+cutname[j]+":alph vs jetE",2);
+ dpsE2[i][j] = dpsf.create(rplnames.get(i)+":"+cutname[j]+":dEoE vs cmE",rplnames.get(i)+":"+cutname[j]+":res vs cmE",2);
+ dpsE3[i][j] = dpsf.create(rplnames.get(i)+":"+cutname[j]+":dEoE vs jetE",rplnames.get(i)+":"+cutname[j]+":res vs jetE",2);
+ }
+ }
+ for(int iE=0;iE<nE;iE++)
+ {
+ ZE = Eevt[iE];
+ System.out.println("Results for cmE = "+ZE+": jetE = "+ZE/2.);
+ for(int k=0;k<rplnames.size();k++)
+ {
+ String rpname = rplnames.get(k);
+ System.out.println(" ReconstructedParticles = "+rpname);
+ String rname = "cmE="+Eevt[iE]+"/"+rpname;
+ System.out.println(" ctcut entries mean rms m90 rms90 alpha90 jetdE/E");
+ for(int i=0;i<cutname.length;i++)
+ {
+ ICloud1D thisc = aida.cloud1D(rname+"/"+cutname[i]+"/dE");
+ double r90 = calc.calculateRMS90(thisc);
+ double m90 = calc.getMEAN90();
+ int ent = thisc.entries();
+ double m = thisc.mean();
+ double r = thisc.rms();
+ double alpha = r90/Math.sqrt(ZE+m90);
+ double dEoE = r90/(ZE+m90);
+ double err = (1.1*r90/Math.sqrt(1.8*ent))/Math.sqrt(ZE+m90);
+ double err2 = (1.1*r90/Math.sqrt(1.8*ent))/(ZE+m90);
+ double jde = r90*Math.sqrt(2.)/(ZE+m90);
+ System.out.println(cutname[i]+" "+ent+" "+m+" "+r+" "+m90+" "+r90+" "+alpha+" "+jde);
+ dpsE[k][i].addPoint();
+ dpsE2[k][i].addPoint();
+ dpsE3[k][i].addPoint();
+ IDataPoint dp = dpsE[k][i].point(iE);
+ IDataPoint dp2 = dpsE2[k][i].point(iE);
+ IDataPoint dp3 = dpsE3[k][i].point(iE);
+ dp.coordinate(0).setValue(ZE/2.);
+ dp.coordinate(1).setValue(alpha);
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ dp2.coordinate(0).setValue(ZE);
+ dp2.coordinate(1).setValue(dEoE);
+ dp2.coordinate(1).setErrorPlus(err2);
+ dp2.coordinate(1).setErrorMinus(err2);
+ dp3.coordinate(0).setValue(ZE/2.);
+ dp3.coordinate(1).setValue(dEoE*Math.sqrt(2.));
+ dp3.coordinate(1).setErrorPlus(err2*Math.sqrt(2.));
+ dp3.coordinate(1).setErrorMinus(err2*Math.sqrt(2.));
+ }
+ IDataPointSet dps = dpsf.create(rname+":alpha90 vs ct",rpname+":alph vs ct",2);
+ IDataPointSet dps2 = dpsf.create(rname+":dEoE vs ct",rpname+":res vs ct",2);
+ IDataPointSet dps5 = dpsf.create(rname+":mean dE vs ct",rpname+":m90 vs ct",2);
+ double ll=0.;
+ for(int i=0;i<ctl.length;i++)
+ {
+ ICloud1D thisc = aida.cloud1D(rname+"/"+ctname[i]+"/dE");
+ double r90 = calc.calculateRMS90(thisc);
+ double m90 = calc.getMEAN90();
+ int ent = thisc.entries();
+ double alpha = r90/Math.sqrt(ZE+m90);
+ double dEoE = r90/(ZE+m90);
+ if(i > 0)ll = ctl[i-1];
+ double x = (ctl[i]+ll)/2.;
+ double err = (1.1*r90/Math.sqrt(1.8*ent))/Math.sqrt(ZE+m90);
+ double err2 = (1.1*r90/Math.sqrt(1.8*ent))/(ZE+m90);
+ double err5 = (1.1*r90/Math.sqrt(.9*ent));
+ dps.addPoint();
+ dps2.addPoint();
+ dps5.addPoint();
+ IDataPoint dp = dps.point(i);
+ IDataPoint dp2 = dps2.point(i);
+ IDataPoint dp5 = dps5.point(i);
+ dp.coordinate(0).setValue(x);
+ dp.coordinate(1).setValue(alpha);
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ dp2.coordinate(0).setValue(x);
+ dp2.coordinate(1).setValue(dEoE);
+ dp2.coordinate(1).setErrorPlus(err2);
+ dp2.coordinate(1).setErrorMinus(err2);
+ dp5.coordinate(0).setValue(x);
+ dp5.coordinate(1).setValue(m90);
+ dp5.coordinate(1).setErrorPlus(err5);
+ dp5.coordinate(1).setErrorMinus(err5);
+ }
+ IDataPointSet dps3 = dpsf.create(rname+":alpha90 vs csct",rname+":alph vs csct",2);
+ IDataPointSet dps4 = dpsf.create(rname+":dEoE vs csct",rname+":res vs csct",2);
+ IDataPointSet dps6 = dpsf.create(rname+":mean dE vs csct",rname+":m90 vs csct",2);
+ ll=1.;
+ for(int i=0;i<cscl.length;i++)
+ {
+ ICloud1D thisc = aida.cloud1D(rname+"/"+cscname[i]+"/dE");
+ double r90 = calc.calculateRMS90(thisc);
+ double m90 = calc.getMEAN90();
+ int ent = thisc.entries();
+ double alpha = r90/Math.sqrt(ZE+m90);
+ double dEoE = r90/(ZE+m90);
+ if(i > 0)ll = cscl[i-1];
+ double x = (cscl[i]+ll)/2.;
+ double err = (1.1*r90/Math.sqrt(1.8*ent))/Math.sqrt(ZE+m90);
+ double err2 = (1.1*r90/Math.sqrt(1.8*ent))/(ZE+m90);
+ double err6 = (1.1*r90/Math.sqrt(.9*ent));
+ dps3.addPoint();
+ dps4.addPoint();
+ dps6.addPoint();
+ IDataPoint dp = dps3.point(i);
+ IDataPoint dp2 = dps4.point(i);
+ IDataPoint dp6 = dps6.point(i);
+ dp.coordinate(0).setValue(x);
+ dp.coordinate(1).setValue(alpha);
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ dp2.coordinate(0).setValue(x);
+ dp2.coordinate(1).setValue(dEoE);
+ dp2.coordinate(1).setErrorPlus(err2);
+ dp2.coordinate(1).setErrorMinus(err2);
+ dp6.coordinate(0).setValue(x);
+ dp6.coordinate(1).setValue(m90);
+ dp6.coordinate(1).setErrorPlus(err6);
+ dp6.coordinate(1).setErrorMinus(err6);
+ }
+ }
+ }
+ super.suspend();
+ }
+}