lcsim/src/org/lcsim/recon/analysis
diff -N RMS90Calculator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RMS90Calculator.java 6 May 2010 17:21:12 -0000 1.1
@@ -0,0 +1,51 @@
+package org.lcsim.recon.analysis;
+import hep.aida.*;
+import java.util.*;
+/*
+ * Calculate RMS90 and Mean90 for a cloud
+ *
+ * @author cassell
+ */
+
+public class RMS90Calculator
+{
+ double rms90;
+ double mean90;
+ public RMS90Calculator()
+ {
+ rms90 = 0.;
+ mean90 = 0.;
+ }
+ public double calculateRMS90(ICloud1D cloud)
+ {
+ int entries = cloud.entries();
+ List<Double> elist = new ArrayList<Double>();
+ for(int i=0;i<entries;i++)
+ {
+ elist.add(cloud.value(i));
+ }
+ rms90 = 999999.;
+ int ntail = (int) (.1*entries);
+ int ncore = entries - ntail;
+ Collections.sort(elist);
+ for(int k=0;k<ntail;k++)
+ {
+ double sv = 0.;
+ double sv2 = 0.;
+ for(int i=k;i<k+ncore;i++)
+ {
+ sv += elist.get(i);
+ sv2 += elist.get(i)*elist.get(i);
+ }
+ double mean = sv/ncore;
+ double rms = Math.sqrt(sv2/ncore - mean*mean);
+ if(rms < rms90)
+ {
+ rms90 = rms;
+ mean90 = mean;
+ }
+ }
+ return rms90;
+ }
+ public double getMEAN90(){return mean90;}
+}
lcsim/src/org/lcsim/recon/analysis
diff -N ZZqqnunuPlots.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ZZqqnunuPlots.java 6 May 2010 17:21:12 -0000 1.1
@@ -0,0 +1,167 @@
+package org.lcsim.recon.analysis;
+
+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 java.util.*;
+import java.io.IOException;
+/*
+ * Make energy and mass plots from a list of ReconstructedParticle for
+ * ZZ -> qqnunu events.
+ *
+ * @author cassell
+ */
+
+public class ZZqqnunuPlots extends Driver
+{
+
+ private AIDA aida = AIDA.defaultInstance();
+ String reconName = "ReconstructedParticles";
+ String pre = "ZZtoqq/";
+ String aidaFN = "ZZevents.aida";
+ int ievt;
+ int nmax = 1000000;
+ double pimass = .1395679;
+ RMS90Calculator calc;
+ double ctcut = 0.97;
+ double ctbarrel = 0.8;
+ public void setCtcut(double c){ctcut = c;}
+ public void setCtbarrel(double c){ctbarrel = c;}
+ public void setAidaFN(String c){aidaFN = c;}
+
+ public ZZqqnunuPlots()
+ {
+ ievt = 0;
+ calc = new RMS90Calculator();
+ }
+
+ protected void process(EventHeader event)
+ {
+ // Find the Z that decays to qq
+ List<MCParticle> mcl = event.get(MCParticle.class, "MCParticle");
+ MCParticle Ztoqq = null;
+ for (MCParticle p : mcl)
+ {
+ if (p.getPDGID() == 23)
+ {
+ if (p.getDaughters().size() == 2)
+ {
+ int id = Math.abs(p.getDaughters().get(0).getPDGID());
+ if ((id == 1) || (id == 2) || (id == 3))
+ {
+ Ztoqq = p;
+ }
+ }
+ }
+ }
+ if (Ztoqq == null)
+ {
+ System.out.println("Could not find Z->qq in event " + ievt);
+ ievt++;
+ return;
+ }
+ double ct0 = Math.abs(Ztoqq.getDaughters().get(0).getMomentum().z())/
+ Ztoqq.getDaughters().get(0).getMomentum().magnitude();
+ double ct1 = Math.abs(Ztoqq.getDaughters().get(1).getMomentum().z())/
+ Ztoqq.getDaughters().get(1).getMomentum().magnitude();
+ super.process(event);
+ //Calculate event energy and mass
+ List<ReconstructedParticle> rpl = event.get(ReconstructedParticle.class, reconName);
+ double Etot = 0.;
+ double px = 0.;
+ double py = 0.;
+ double pz = 0.;
+ for (ReconstructedParticle p : rpl)
+ {
+ // Put in fix for earlier reconstruction error
+ double pE = p.getEnergy();
+ Hep3Vector mom = p.getMomentum();
+ if (Double.isNaN(pE))
+ {
+ if (p.getCharge() == 0.)
+ {
+ System.out.println("Skipping event " + ievt + " with undefined neutral Energy");
+ ievt++;
+ return;
+ } else
+ {
+ mom = new BasicHep3Vector(p.getTracks().get(0).getMomentum());
+ pE = Math.sqrt(mom.magnitudeSquared() + pimass * pimass);
+ }
+ }
+ px += mom.x();
+ py += mom.y();
+ pz += mom.z();
+ Etot += pE;
+ }
+ double M = Math.sqrt(Etot * Etot - px * px - py * py - pz * pz);
+ double dE = Etot - Ztoqq.getEnergy();
+ double dM = M - Ztoqq.getMass();
+ String pre2 = "All events/";
+ aida.cloud1D(pre+pre2+"Gen event Energy",nmax).fill(Ztoqq.getEnergy());
+ aida.cloud1D(pre+pre2+"Gen event Mass",nmax).fill(Ztoqq.getMass());
+ aida.cloud1D(pre+pre2+"Recon event Energy",nmax).fill(Etot);
+ aida.cloud1D(pre+pre2+"Recon event Mass",nmax).fill(M);
+ aida.cloud1D(pre+pre2+"Delta Energy",nmax).fill(dE);
+ aida.cloud1D(pre+pre2+"Delta Mass",nmax).fill(dM);
+ if( (ct0 < ctcut)&&(ct1 < ctcut) )
+ {
+ pre2 = "Events ct < "+ctcut+"/";
+ aida.cloud1D(pre+pre2+"Gen event Energy",nmax).fill(Ztoqq.getEnergy());
+ aida.cloud1D(pre+pre2+"Gen event Mass",nmax).fill(Ztoqq.getMass());
+ aida.cloud1D(pre+pre2+"Recon event Energy",nmax).fill(Etot);
+ aida.cloud1D(pre+pre2+"Recon event Mass",nmax).fill(M);
+ aida.cloud1D(pre+pre2+"Delta Energy",nmax).fill(dE);
+ aida.cloud1D(pre+pre2+"Delta Mass",nmax).fill(dM);
+ pre2 = "Split BE/";
+ if(ct0 < ctbarrel)
+ {
+ if(ct1 < ctbarrel)pre2 = "Both Barrel/";
+ }
+ else if(ct1 >= ctbarrel)pre2 = "Both Endcap/";
+ aida.cloud1D(pre+pre2+"Gen event Energy",nmax).fill(Ztoqq.getEnergy());
+ aida.cloud1D(pre+pre2+"Gen event Mass",nmax).fill(Ztoqq.getMass());
+ aida.cloud1D(pre+pre2+"Recon event Energy",nmax).fill(Etot);
+ aida.cloud1D(pre+pre2+"Recon event Mass",nmax).fill(M);
+ aida.cloud1D(pre+pre2+"Delta Energy",nmax).fill(dE);
+ aida.cloud1D(pre+pre2+"Delta Mass",nmax).fill(dM);
+ }
+ else
+ {
+ pre2 = "Events ct >= "+ctcut+"/";
+ aida.cloud1D(pre+pre2+"Gen event Energy",nmax).fill(Ztoqq.getEnergy());
+ aida.cloud1D(pre+pre2+"Gen event Mass",nmax).fill(Ztoqq.getMass());
+ aida.cloud1D(pre+pre2+"Recon event Energy",nmax).fill(Etot);
+ aida.cloud1D(pre+pre2+"Recon event Mass",nmax).fill(M);
+ aida.cloud1D(pre+pre2+"Delta Energy",nmax).fill(dE);
+ aida.cloud1D(pre+pre2+"Delta Mass",nmax).fill(dM);
+ }
+ ievt++;
+ }
+ protected void suspend()
+ {
+ try{
+ AIDA.defaultInstance().saveAs(aidaFN);
+ }
+ catch(IOException ioe1)
+ {
+ ioe1.printStackTrace();
+ }
+ String[] fold = {"All events/","Events ct < "+ctcut+"/","Events ct >= "+ctcut+"/",
+ "Split BE/","Both Barrel/","Both Endcap/"};
+ String[] name = {"Delta Energy","Delta Mass"};
+ for(int i=0;i<fold.length;i++)
+ {
+ for(int j=0;j<name.length;j++)
+ {
+ String pp = pre+fold[i]+name[j];
+ double rms90 = calc.calculateRMS90(aida.cloud1D(pp));
+ double mean90 = calc.getMEAN90();
+ System.out.println(pp+": mean90 = "+mean90+", rms90 = "+rms90);
+ }
+ }
+ }
+}