lcsim/src/org/lcsim/contrib/SteveMagill
diff -N PFASelect.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PFASelect.java 18 May 2007 19:59:11 -0000 1.1
@@ -0,0 +1,170 @@
+package org.lcsim.contrib.SteveMagill;
+
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.util.Driver;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.jet.*;
+
+/*
+ * PFASelect includes fiducial cut to select region for PFA analysis
+ */
+
+public class PFASelect extends Driver
+{
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ public PFASelect()
+ {
+ add(new PFATemplate());
+
+ // now find jets with real Reco Particles
+ JetDriver PFAjet = new JetDriver();
+ PFAjet.setInputCollectionName("AllRecoParticles");
+ PFAjet.setOutputCollectionName("PFAJets");
+ JetFinder tworealjets = new FixNumberOfJetsFinder(2);
+ PFAjet.setFinder(tworealjets);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder rktjet = new DurhamJetFinder(0.0005);
+// PFAjet.setFinder(rktjet);
+ add(PFAjet);
+
+ // make perfect jets out of perfect particles
+ JetDriver j = new JetDriver();
+ j.setInputCollectionName("PerfectRecoParticles");
+ j.setOutputCollectionName("PerfectJets");
+ JetFinder twojet = new FixNumberOfJetsFinder(2);
+ j.setFinder(twojet);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder ktjet = new DurhamJetFinder(0.0005);
+// j.setFinder(ktjet);
+ add(j);
+
+ // make plots etc. from results
+ PPFAPlotDriver ppfaplots = new PPFAPlotDriver();
+ ppfaplots.setMCPartName("CheatReconstructedParticles");
+ ppfaplots.setPerfJetName("PerfectJets");
+ ppfaplots.setPFAJetName("PFAJets");
+ ppfaplots.setPerfPartName("PerfectRecoParticles");
+ ppfaplots.setPFAPartName("AllRecoParticles");
+ add(ppfaplots);
+ }
+
+ protected void process(EventHeader event)
+ {
+ // Only process events with both jets costheta < costhetaMAX
+
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+ boolean keepit = true;
+ double Zmass = 0.;
+ int ievt = 0;
+ double costmax = 0.8;
+ 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();
+// System.out.println("Z Mass " + Zmass);
+ Hep3Vector pp = p.getMomentum();
+ double cost = Math.abs(pp.z()/pp.magnitude());
+ if(cost > costmax)
+ {
+ keepit = false;
+ continue;
+ }
+ }
+ }
+ }
+
+ if(!keepit)
+ {
+ ievt++;
+ return;
+ }
+ super.process(event);
+ // add some plots here - get jet stuff
+ List<ReconstructedParticle> perfjets = event.get(ReconstructedParticle.class,"PerfectJets");
+ List<ReconstructedParticle> pfajets = event.get(ReconstructedParticle.class,"PFAJets");
+
+ // here's some jet stuff for perfect jets
+ double jESum = 0;
+ double jPXSum = 0;
+ double jPYSum = 0;
+ double jPZSum = 0;
+ double djm = 0;
+ double[] jpz = new double[2];
+ double[] jp = new double[2];
+ int j = 0;
+ aida.cloud1D("Number of Perfect Jets").fill(perfjets.size());
+ if (perfjets.size() == 2)
+ {
+ // calculate dijet mass - don't we have something that does this?
+ for (ReconstructedParticle perfjet : perfjets)
+ {
+ aida.cloud1D("Perfect Jet Energy").fill(perfjet.getEnergy());
+ aida.cloud1D("Perfect Jet Mass").fill(perfjet.getMass());
+// jESum += perfjet.getEnergy();
+ Hep3Vector j3v = perfjet.getMomentum();
+ jPXSum += j3v.x();
+ jPYSum += j3v.y();
+ jPZSum += j3v.z();
+ jp[j] = Math.sqrt(j3v.x()*j3v.x()+j3v.y()*j3v.y()+j3v.z()*j3v.z());
+ // force jet mass 0
+ jESum += jp[j];
+ aida.cloud1D("Perfect Jet Energy zero mass").fill(jp[j]);
+ jpz[j] = j3v.z();
+ j++;
+ }
+ double jE2 = jESum*jESum;
+ double jP2 = jPXSum*jPXSum+jPYSum*jPYSum+jPZSum*jPZSum;
+ djm += Math.sqrt(jE2-jP2);
+ aida.cloud1D("Perfect Dijet Mass").fill(djm);
+ aida.cloud1D("Diff Perfect Dijet Mass").fill((djm-Zmass)/Math.sqrt(Zmass));
+ }
+
+ // now real PFA jets
+ double rjESum = 0;
+ double rjPXSum = 0;
+ double rjPYSum = 0;
+ double rjPZSum = 0;
+ double rdjm = 0;
+ double[] jrpz = new double[2];
+ double[] jrp = new double[2];
+ int k = 0;
+ aida.cloud1D("Number of PFA Jets").fill(pfajets.size());
+ if (pfajets.size() == 2)
+ {
+ // calculate dijet mass - don't we have something that does this?
+ for (ReconstructedParticle pfajet : pfajets)
+ {
+ aida.cloud1D("PFA Jet Energy").fill(pfajet.getEnergy());
+ aida.cloud1D("PFA Jet Mass").fill(pfajet.getMass());
+// rjESum += pfajet.getEnergy();
+ Hep3Vector rj3v = pfajet.getMomentum();
+ rjPXSum += rj3v.x();
+ rjPYSum += rj3v.y();
+ rjPZSum += rj3v.z();
+ jrp[k] = Math.sqrt(rj3v.x()*rj3v.x()+rj3v.y()*rj3v.y()+rj3v.z()*rj3v.z());
+ // force jet mass 0
+ rjESum += jrp[k];
+ aida.cloud1D("PFA Jet Energy zero mass").fill(jrp[k]);
+ jrpz[k] = rj3v.z();
+ k++;
+ }
+ double rjE2 = rjESum*rjESum;
+ double rjP2 = rjPXSum*rjPXSum+rjPYSum*rjPYSum+rjPZSum*rjPZSum;
+ rdjm += Math.sqrt(rjE2-rjP2);
+ aida.cloud1D("PFA Dijet Mass").fill(rdjm);
+ aida.cloud1D("Diff PFA Dijet Mass").fill((rdjm-Zmass)/Math.sqrt(Zmass));
+ aida.cloud1D("Difference PFAJet and PerfJet DiJet Mass").fill(rdjm-djm);
+ }
+ }
+}