lcsim/src/org/lcsim/contrib/SteveMagill
diff -N PFASelect.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PFASelect.java 20 Aug 2008 20:18:31 -0000 1.3
@@ -0,0 +1,204 @@
+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 HighPTemplate());
+// add(new MipRPhoTSCNeuTemplate());
+// add(new MHPTrShPNTemplate());
+// 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);
+// JetFinder sixrealjets = new FixNumberOfJetsFinder(6);
+// PFAjet.setFinder(sixrealjets);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder rktjet = new DurhamJetFinder(0.0003);
+// 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);
+// JetFinder sixjet = new FixNumberOfJetsFinder(6);
+// j.setFinder(sixjet);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder ktjet = new DurhamJetFinder(0.0003);
+// 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.;
+ double qqbarE = 0.;
+ int ievt = 0;
+ double costmax = 0.8; // 0.8 is barrel
+ 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();
+ qqbarE = p.getParents().get(0).getEnergy();
+// 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");
+ List<ReconstructedParticle> rpList = event.get(ReconstructedParticle.class,"AllRecoParticles");
+
+ // check RP stuff
+ double rpESum = 0.;
+ double rpPxSum = 0.;
+ double rpPySum = 0.;
+ double rpPzSum = 0.;
+ for(ReconstructedParticle rp : rpList)
+ {
+ aida.cloud1D("AllRecoParticle Energy").fill(rp.getEnergy());
+ rpESum += rp.getEnergy();
+ Hep3Vector rpmom = rp.getMomentum();
+ rpPxSum += rpmom.x();
+ rpPySum += rpmom.y();
+ rpPzSum += rpmom.z();
+ }
+ double rpPSum = Math.sqrt(rpPxSum*rpPxSum+rpPySum*rpPySum+rpPzSum*rpPzSum);
+ aida.cloud1D("AllRecoParticle Diff ESum").fill(rpESum-qqbarE);
+ aida.cloud1D("AllRecoParticle PSum").fill(rpPSum);
+ aida.cloud1D("AllRecoParticle Diff Mass").fill(Math.sqrt(rpESum*rpESum-rpPSum*rpPSum)-Zmass);
+
+ // 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());
+ aida.cloud1D("Number of particles in Jet").fill(perfjet.getParticles().size());
+ 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("Z Mass").fill(Zmass);
+ aida.cloud1D("Perfect Dijet Mass").fill(djm);
+ aida.cloud1D("Perfect Dijet Energy").fill(Math.sqrt(jE2));
+ aida.cloud1D("Diff Perfect Dijet Mass").fill(djm-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("PFA Dijet Energy").fill(Math.sqrt(rjE2));
+ aida.cloud1D("Diff PFA Dijet Mass").fill(rdjm-Zmass);
+ aida.cloud1D("Difference PFAJet and PerfJet DiJet Mass").fill(rdjm-djm);
+ }
+ }
+}