lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat
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 ++;
+ }
+}