lcsim/src/org/lcsim/recon/analysis

diff -u -r1.1 -r1.2
--- ZZqqnunuPlots.java	6 May 2010 17:21:12 -0000	1.1
+++ ZZqqnunuPlots.java	18 Mar 2011 12:30:18 -0000	1.2
@@ -3,6 +3,7 @@

 import org.lcsim.util.Driver;
 import org.lcsim.event.EventHeader;
 import org.lcsim.event.MCParticle;

+import org.lcsim.event.Track;

 import org.lcsim.event.ReconstructedParticle;
 import hep.physics.vec.*;
 import org.lcsim.util.aida.AIDA;

@@ -31,6 +32,7 @@

     public void setCtcut(double c){ctcut = c;}
     public void setCtbarrel(double c){ctbarrel = c;}
     public void setAidaFN(String c){aidaFN = c;}

+    public void setReconName(String c){reconName = c;}

 
     public ZZqqnunuPlots()
     {

@@ -42,6 +44,7 @@

     {
         // Find the Z that decays to qq
         List<MCParticle> mcl = event.get(MCParticle.class, "MCParticle");

+        List<Track> tl = new ArrayList<Track>(event.get(Track.class, "Tracks"));

         MCParticle Ztoqq = null;
         for (MCParticle p : mcl)
         {

@@ -74,6 +77,13 @@

         double px = 0.;
         double py = 0.;
         double pz = 0.;

+        double Ech = 0.;
+        double Enh = 0.;
+        double Eph = 0.;
+        double Echcl = 0.;
+        int nch = 0;
+        int nnh = 0;
+        int nph = 0;

         for (ReconstructedParticle p : rpl)
         {
             // Put in fix for earlier reconstruction error

@@ -96,10 +106,41 @@

             py += mom.y();
             pz += mom.z();
             Etot += pE;

+            if(p.getCharge() == 0.)
+            {
+               if(p.getMass() == 0.)
+               {
+                  Eph += pE;
+                  nph++;
+               }
+               else
+               {
+                  Enh += pE;
+                  nnh++;
+               }
+            }
+            else
+            {
+               Ech += pE;
+               nch++;
+               tl.remove(p.getTracks().get(0));
+               double clE = 0.;
+               if(p.getClusters().size() > 0)clE = p.getClusters().get(0).getEnergy();
+               Echcl += clE;
+            }
+        }
+        int nuut = tl.size();
+        double uutE = 0.;
+        for(Track t:tl)
+        {
+           Hep3Vector mom = new BasicHep3Vector(t.getMomentum());
+           uutE += Math.sqrt(mom.magnitudeSquared()+pimass*pimass);

         }
         double M = Math.sqrt(Etot * Etot - px * px - py * py - pz * pz);
         double dE = Etot - Ztoqq.getEnergy();
         double dM = M - Ztoqq.getMass();

+        double dEoE = dE/Ztoqq.getEnergy();
+        double dMoM = dM/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());

@@ -107,6 +148,17 @@

         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);

+        aida.cloud1D(pre+pre2+"DE over E",nmax).fill(dEoE);
+        aida.cloud1D(pre+pre2+"DM over M",nmax).fill(dMoM);
+        aida.cloud1D(pre+pre2+"Recon charged E per event",nmax).fill(Ech);
+        aida.cloud1D(pre+pre2+"Recon nhadron E per event",nmax).fill(Enh);
+        aida.cloud1D(pre+pre2+"Recon photon E per event",nmax).fill(Eph);
+        aida.cloud1D(pre+pre2+"# charged per event",nmax).fill(nch);
+        aida.cloud1D(pre+pre2+"# nhadron per event",nmax).fill(nnh);
+        aida.cloud1D(pre+pre2+"# photon per event",nmax).fill(nph);
+        aida.cloud1D(pre+pre2+"# unused tracks per event",nmax).fill(nuut);
+        aida.cloud1D(pre+pre2+"unused track E per event",nmax).fill(uutE);
+        aida.cloud1D(pre+pre2+"Recon charged cluster E per event",nmax).fill(Echcl);

         if( (ct0 < ctcut)&&(ct1 < ctcut) )
         {
             pre2 = "Events ct < "+ctcut+"/";

@@ -116,6 +168,17 @@

             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);

+            aida.cloud1D(pre+pre2+"DE over E",nmax).fill(dEoE);
+            aida.cloud1D(pre+pre2+"DM over M",nmax).fill(dMoM);
+            aida.cloud1D(pre+pre2+"Recon charged E per event",nmax).fill(Ech);
+            aida.cloud1D(pre+pre2+"Recon nhadron E per event",nmax).fill(Enh);
+            aida.cloud1D(pre+pre2+"Recon photon E per event",nmax).fill(Eph);
+            aida.cloud1D(pre+pre2+"# charged per event",nmax).fill(nch);
+            aida.cloud1D(pre+pre2+"# nhadron per event",nmax).fill(nnh);
+            aida.cloud1D(pre+pre2+"# photon per event",nmax).fill(nph);
+            aida.cloud1D(pre+pre2+"# unused tracks per event",nmax).fill(nuut);
+            aida.cloud1D(pre+pre2+"unused track E per event",nmax).fill(uutE);
+            aida.cloud1D(pre+pre2+"Recon charged cluster E per event",nmax).fill(Echcl);

             pre2 = "Split BE/";
             if(ct0 < ctbarrel)
             {

@@ -128,6 +191,36 @@

             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);

+            aida.cloud1D(pre+pre2+"DE over E",nmax).fill(dEoE);
+            aida.cloud1D(pre+pre2+"DM over M",nmax).fill(dMoM);
+            aida.cloud1D(pre+pre2+"Recon charged E per event",nmax).fill(Ech);
+            aida.cloud1D(pre+pre2+"Recon nhadron E per event",nmax).fill(Enh);
+            aida.cloud1D(pre+pre2+"Recon photon E per event",nmax).fill(Eph);
+            aida.cloud1D(pre+pre2+"# charged per event",nmax).fill(nch);
+            aida.cloud1D(pre+pre2+"# nhadron per event",nmax).fill(nnh);
+            aida.cloud1D(pre+pre2+"# photon per event",nmax).fill(nph);
+            aida.cloud1D(pre+pre2+"# unused tracks per event",nmax).fill(nuut);
+            aida.cloud1D(pre+pre2+"unused track E per event",nmax).fill(uutE);
+            aida.cloud1D(pre+pre2+"Recon charged cluster E per event",nmax).fill(Echcl);
+            String pre3 = "+/";
+            if(Ztoqq.getDaughters().get(0).getMomentum().z() < 0.)pre3 = "-/";
+            aida.cloud1D(pre+pre2+pre3+"Gen event Energy",nmax).fill(Ztoqq.getEnergy());
+            aida.cloud1D(pre+pre2+pre3+"Gen event Mass",nmax).fill(Ztoqq.getMass());
+            aida.cloud1D(pre+pre2+pre3+"Recon event Energy",nmax).fill(Etot);
+            aida.cloud1D(pre+pre2+pre3+"Recon event Mass",nmax).fill(M);
+            aida.cloud1D(pre+pre2+pre3+"Delta Energy",nmax).fill(dE);
+            aida.cloud1D(pre+pre2+pre3+"Delta Mass",nmax).fill(dM);
+            aida.cloud1D(pre+pre2+pre3+"DE over E",nmax).fill(dEoE);
+            aida.cloud1D(pre+pre2+pre3+"DM over M",nmax).fill(dMoM);
+            aida.cloud1D(pre+pre2+pre3+"Recon charged E per event",nmax).fill(Ech);
+            aida.cloud1D(pre+pre2+pre3+"Recon nhadron E per event",nmax).fill(Enh);
+            aida.cloud1D(pre+pre2+pre3+"Recon photon E per event",nmax).fill(Eph);
+            aida.cloud1D(pre+pre2+pre3+"# charged per event",nmax).fill(nch);
+            aida.cloud1D(pre+pre2+pre3+"# nhadron per event",nmax).fill(nnh);
+            aida.cloud1D(pre+pre2+pre3+"# photon per event",nmax).fill(nph);
+            aida.cloud1D(pre+pre2+pre3+"# unused tracks per event",nmax).fill(nuut);
+            aida.cloud1D(pre+pre2+pre3+"unused track E per event",nmax).fill(uutE);
+            aida.cloud1D(pre+pre2+pre3+"Recon charged cluster E per event",nmax).fill(Echcl);

         }
         else
         {

@@ -138,6 +231,17 @@

             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);

+            aida.cloud1D(pre+pre2+"DE over E",nmax).fill(dEoE);
+            aida.cloud1D(pre+pre2+"DM over M",nmax).fill(dMoM);
+            aida.cloud1D(pre+pre2+"Recon charged E per event",nmax).fill(Ech);
+            aida.cloud1D(pre+pre2+"Recon nhadron E per event",nmax).fill(Enh);
+            aida.cloud1D(pre+pre2+"Recon photon E per event",nmax).fill(Eph);
+            aida.cloud1D(pre+pre2+"# charged per event",nmax).fill(nch);
+            aida.cloud1D(pre+pre2+"# nhadron per event",nmax).fill(nnh);
+            aida.cloud1D(pre+pre2+"# photon per event",nmax).fill(nph);
+            aida.cloud1D(pre+pre2+"# unused tracks per event",nmax).fill(nuut);
+            aida.cloud1D(pre+pre2+"unused track E per event",nmax).fill(uutE);
+            aida.cloud1D(pre+pre2+"Recon charged cluster E per event",nmax).fill(Echcl);

         }
         ievt++;
     }

@@ -161,6 +265,17 @@

                 double rms90 = calc.calculateRMS90(aida.cloud1D(pp));
                 double mean90 = calc.getMEAN90();
                 System.out.println(pp+": mean90 = "+mean90+", rms90 = "+rms90);

+                if(i>2)
+                {
+                   pp = pre+fold[i]+"+/"+name[j];
+                   rms90 = calc.calculateRMS90(aida.cloud1D(pp));
+                   mean90 = calc.getMEAN90();
+                   System.out.println(pp+": mean90 = "+mean90+", rms90 = "+rms90);
+                   pp = pre+fold[i]+"-/"+name[j];
+                   rms90 = calc.calculateRMS90(aida.cloud1D(pp));
+                   mean90 = calc.getMEAN90();
+                   System.out.println(pp+": mean90 = "+mean90+", rms90 = "+rms90);
+                }

             }
         }
     }