lcsim/src/org/lcsim/contrib/SteveMagill
diff -N PPFAPlotDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PPFAPlotDriver.java 23 Apr 2007 20:05:34 -0000 1.1
@@ -0,0 +1,219 @@
+package org.lcsim.contrib.SteveMagill;
+
+import java.util.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.geometry.util.CalorimeterIDDecoder;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.event.Track;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.event.ParticleID;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
+
+/**
+ * Defines perfect particle flow objects - perfect charged particles, perfect photons and perfect neutral hadrons
+ */
+
+public class PPFAPlotDriver extends Driver
+{
+ public PPFAPlotDriver() {
+ }
+
+ private String _perfname;
+ private String _mcname;
+ private String _perfjetname;
+ private String _pfajetname;
+ private AIDA aida = AIDA.defaultInstance();
+ private boolean perfPFAD = true;
+
+ public void process(EventHeader event)
+ {
+ List<ReconstructedParticle> perfrps = event.get(ReconstructedParticle.class,_perfname);
+ List<MCReconstructedParticle> mcrps = event.get(MCReconstructedParticle.class,_mcname);
+ List<ReconstructedParticle> perfjets = event.get(ReconstructedParticle.class,_perfjetname);
+ List<ReconstructedParticle> pfajets = event.get(ReconstructedParticle.class, _pfajetname);
+ List<CheatCluster> chclusters = event.get(CheatCluster.class,"PerfectCheatClusters");
+ List<BasicCluster> phobcls = event.get(BasicCluster.class,"PhotonBClusters");
+ List<BasicCluster> phoeccls = event.get(BasicCluster.class,"PhotonECClusters");
+ List<ReconstructedParticle> ppjets = event.get(ReconstructedParticle.class,"PerfectJets");
+
+// List<CheatCluster> chHADclusters = event.get(CheatCluster.class,"HcalBarrDigiHitsCheatClusters");
+
+ // 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;
+ if (perfPFAD) 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)
+ {
+ if (perfPFAD) aida.cloud1D("Perfect Jet Energy").fill(perfjet.getEnergy());
+ 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());
+ jpz[j] = j3v.z();
+ j++;
+ }
+ double jE2 = jESum*jESum;
+ double jP2 = jPXSum*jPXSum+jPYSum*jPYSum+jPZSum*jPZSum;
+ djm += Math.sqrt(jE2-jP2);
+ if (perfPFAD && jpz[0]/jp[0]<0.8 && jpz[1]/jp[1]<0.8) aida.cloud1D("Perf Barrel DiJet Mass").fill(djm);
+ if (perfPFAD) aida.cloud1D("Perfect Dijet Mass").fill(djm);
+ if (perfPFAD) aida.cloud1D("Diff Perfect Dijet Mass").fill((djm-91)/9.5);
+ }
+
+ // 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;
+ if (perfPFAD) 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)
+ {
+ if (perfPFAD) aida.cloud1D("PFA Jet Energy").fill(pfajet.getEnergy());
+ 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());
+ jrpz[k] = rj3v.z();
+ k++;
+ }
+ double rjE2 = rjESum*rjESum;
+ double rjP2 = rjPXSum*rjPXSum+rjPYSum*rjPYSum+rjPZSum*rjPZSum;
+ rdjm += Math.sqrt(rjE2-rjP2);
+ if (perfPFAD) aida.cloud1D("PFA Dijet Mass").fill(rdjm);
+ if (perfPFAD) aida.cloud1D("Diff PFA Dijet Mass").fill((rdjm-91)/9.5);
+ if (perfPFAD && jrpz[0]/jrp[0]<0.8 && jrpz[1]/jrp[1]<0.8) aida.cloud1D("PFA Barrel DiJet Mass").fill(rdjm);
+ if (perfPFAD && rdjm>0 && djm>0) aida.cloud1D("Difference PFAJet and PerfJet DiJet Mass").fill((rdjm-djm)/Math.sqrt(djm));
+ }
+
+ // Check some particle comparisons - photons first
+ int NPFAPho = phobcls.size()+phoeccls.size();
+ double PFAPhoClE = 0;
+ for (BasicCluster phobcl : phobcls)
+ {
+ PFAPhoClE += phobcl.getEnergy()*1.015*1.013;
+ }
+ for (BasicCluster phoeccl : phoeccls)
+ {
+ PFAPhoClE += phoeccl.getEnergy()*1.015*1.013;
+ }
+
+ // now some stuff to compare perfect reco particles to true particles
+ double perfpho = 0;
+ double perfchpE = 0;
+ double perfchpP = 0;
+ double perfneu = 0;
+ for (ReconstructedParticle perfrp : perfrps)
+ {
+ double perfpCh = perfrp.getCharge();
+// ParticleID perfpID = perfrp.getParticleIDUsed();
+ int perfpID = Math.abs(perfrp.getParticleIDUsed().getPDG());
+// System.out.println("Particle ID = " + perfpID);
+ double perfpE = perfrp.getEnergy();
+ Hep3Vector perfpP = perfrp.getMomentum();
+ double pP = perfpP.magnitude();
+ if (perfpCh == 0)
+ {
+ if (perfpID == 22)
+ {
+ perfpho += perfpE;
+ }else
+ {
+ perfneu += perfpE;
+ }
+ } else
+ {
+ perfchpE += perfpE;
+ perfchpP += pP;
+ }
+ for (MCReconstructedParticle mcrp : mcrps)
+ {
+ MCParticle mcp = mcrp.getMCParticle();
+ }
+ }
+
+ // Check the cheated clusters also
+ double chclusE = 0;
+ int NPerfPho = 0;
+ double PerfPhoClE = 0;
+ for (CheatCluster chclus : chclusters)
+ {
+ chclusE += chclus.getEnergy();
+ if (chclus.getMCParticle().getCharge() == 0 && chclus.getMCParticle().getPDGID() == 22)
+ {
+ NPerfPho++;
+ PerfPhoClE += chclus.getEnergy();
+ }
+ if (chclus.getMCParticle().getCharge() != 0)
+ {
+ double chclchpE = chclus.getEnergy();
+ double chclchpP = chclus.getMCParticle().getMomentum().magnitude();
+ if (perfPFAD) aida.cloud1D("Cheated Cluster Charge E over P").fill(chclchpE/chclchpP);
+ }
+ }
+ if (perfPFAD) aida.cloud1D("Cheated CLuster E").fill(chclusE);
+ if (perfPFAD) aida.cloud1D("Difference NPFAPhotons NPerfPhotons").fill(NPFAPho-NPerfPho);
+ if (perfPFAD) aida.cloud1D("Difference PFAPhoE PerfPhoE").fill(PFAPhoClE-PerfPhoClE);
+
+ double TPESum = perfpho+perfchpP+perfneu;
+ if (perfPFAD) aida.cloud1D("Total Perfect ESum").fill(TPESum);
+ if (perfPFAD) aida.cloud1D("Perfect Photon ESum").fill(perfpho);
+ if (perfPFAD) aida.cloud1D("Perfect Charge P ESum").fill(perfchpE);
+ if (perfPFAD) aida.cloud1D("Perfect Charge P PSum").fill(perfchpP);
+ if (perfPFAD) aida.cloud1D("Perfect E over P for Charge P").fill(perfchpP/perfchpE);
+ if (perfPFAD) aida.cloud1D("Perfect Neutral H ESum").fill(perfneu);
+ if (perfPFAD) aida.cloud1D("Perfect Charge P + Photon ESum").fill(perfchpP+perfpho);
+ if (perfPFAD) aida.cloud1D("Perf Pho ESum Ratio").fill(perfpho/TPESum);
+ if (perfPFAD) aida.cloud1D("Perf Ch P ESum Ratio").fill(perfchpP/TPESum);
+ if (perfPFAD) aida.cloud1D("Perf Neu H ESum Ratio").fill(perfneu/TPESum);
+ }
+
+ public void setPerfPartName(String perfname)
+ {
+ _perfname = perfname;
+ }
+ public void setMCPartName(String mcname)
+ {
+ _mcname = mcname;
+ }
+ public void setPerfJetName(String perfjetname)
+ {
+ _perfjetname = perfjetname;
+ }
+ public void setPFAJetName(String pfajetname)
+ {
+ _pfajetname = pfajetname;
+ }
+
+}
