Print

Print


Commit in lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat on MAIN
example2JetAnalysis.java+103added 1.1
Example driver to generate pattern recognition cheating reconstructed particles, and then find 2 jets and analyse

lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat
example2JetAnalysis.java added at 1.1
diff -N example2JetAnalysis.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ example2JetAnalysis.java	15 Feb 2007 23:03:42 -0000	1.1
@@ -0,0 +1,103 @@
+import java.util.*;
+import hep.aida.*;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.jet.*;
+import hep.physics.vec.*;
+
+/**
+ * Test analysis using reconstructed particles
+ */
+public class example2JetAnalysis extends Driver
+{
+   private AIDA aida = AIDA.defaultInstance();
+   String reconname = "PPRReconParticles";
+   String jetlistname = "PPRJetReconstructedParticles";
+   int ievt;
+   public example2JetAnalysis()
+   {
+// 
+//  Create a List of perfect pattern recognition reconstructed particles
+//
+      PPRReconDriver rd = new PPRReconDriver();
+      rd.setOutputName(reconname);
+      add(rd);
+//
+//  Ask the jet finder to find 2 jets
+//
+      JetDriver j = new JetDriver();
+      j.setInputCollectionName(reconname);
+      j.setOutputCollectionName(jetlistname);
+      JetFinder twojet = new FixNumberOfJetsFinder(2);
+      j.setFinder(twojet);
+      add(j);
+      
+      ievt = 0;
+   }
+   protected void process(EventHeader event)
+   {
+       super.process(event);
+       List<ReconstructedParticle>jets = event.get(ReconstructedParticle.class,jetlistname);
+       aida.cloud1D("# jets per event").fill(jets.size());
+       double esum = 0.;
+       double pxsum = 0.;
+       double pysum = 0.;
+       double pzsum = 0.;
+       double chargesum = 0.;
+       double[] je = new double[jets.size()];
+       Hep3Vector[] jp = new Hep3Vector[jets.size()];
+       int ind = 0;
+       for(ReconstructedParticle jet:jets)
+       {
+          esum += jet.getEnergy();
+          je[ind] = jet.getEnergy();
+          Hep3Vector v = jet.getMomentum();
+          jp[ind] = v;
+          pxsum += v.x();
+          pysum += v.y();
+          pzsum += v.z();
+          chargesum += jet.getCharge();
+          ind++;
+       }
+       double psum = Math.sqrt(pxsum*pxsum+pysum*pysum+pzsum*pzsum);
+       double evtmass = Math.sqrt(esum*esum - psum*psum);
+       aida.cloud1D(" Event energy").fill(esum);
+       aida.cloud1D("Event momentum").fill(psum);
+       aida.cloud1D("Event charge").fill(chargesum);
+       aida.cloud1D("Event mass").fill(evtmass);
+       for(int i=0;i<jets.size();i++)
+       {
+          double ct = jp[i].z()/Math.sqrt(jp[i].x()*jp[i].x() + jp[i].y()*jp[i].y() + jp[i].z()*jp[i].z());
+          double jm = Math.sqrt(jp[i].x()*jp[i].x() + jp[i].y()*jp[i].y() + jp[i].z()*jp[i].z());
+          aida.cloud1D("Jet energy").fill(je[i]);
+          aida.cloud1D("Jet momentum").fill(jm);
+          aida.cloud1D("Jet mass").fill(Math.sqrt(je[i]*je[i] - jm*jm));
+          aida.cloud1D("cosTheta jet").fill(jp[i].z()/jm);
+          if(Math.abs(ct) < .8)
+          {
+             aida.cloud1D("Cut Jet energy").fill(je[i]);
+             aida.cloud1D("Cut Jet momentum").fill(jm);
+             aida.cloud1D("Cut Jet mass").fill(Math.sqrt(je[i]*je[i] - jm*jm));
+          }
+          for(int j=i+1;j<jets.size();j++)
+          {
+             double ct1 = jp[j].z()/Math.sqrt(jp[j].x()*jp[j].x() + jp[j].y()*jp[j].y() + jp[j].z()*jp[j].z());
+             double masssq = (je[i]+je[j])*(je[i]+je[j]) - (jp[i].x()+jp[j].x())*(jp[i].x()+jp[j].x()) - (jp[i].y()+jp[j].y())*(jp[i].y()+jp[j].y())
+                                               - (jp[i].z()+jp[j].z())*(jp[i].z()+jp[j].z());
+             double sign = 1.;
+             if(masssq < 0.)sign = -1.;
+             aida.cloud1D("Jet - Jet mass").fill(sign*Math.sqrt(sign*masssq));
+             aida.cloud2D("cosTh1 vs cosTh2").fill(ct1,ct);
+             if( (Math.abs(ct) < .8)&&(Math.abs(ct1) < .8) )
+             {
+                aida.cloud1D("Cut Jet - Jet mass").fill(sign*Math.sqrt(sign*masssq));
+             }
+          }
+       }
+       ievt ++;
+   }
+}
CVSspam 0.2.8