lcsim/src/org/lcsim/contrib/SteveMagill
diff -N PhoCompDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PhoCompDriver.java 29 May 2008 21:30:15 -0000 1.1
@@ -0,0 +1,188 @@
+package org.lcsim.contrib.compile.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.spacegeom.*;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.event.ParticleID;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
+import org.lcsim.event.ParticleID;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseParticleID;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticleType;
+
+/**
+ * Defines perfect particle flow objects - perfect charged particles, perfect photons and perfect neutral hadrons
+ */
+
+public class PhoCompDriver extends Driver
+{
+ public PhoCompDriver() {
+ }
+
+ private String _phnames;
+ private AIDA aida = AIDA.defaultInstance();
+ private boolean phoD = true;
+
+ public void process(EventHeader event)
+ {
+ List<CheatCluster> chclusters = event.get(CheatCluster.class,"PerfectCheatClusters");
+ List<BasicCluster> phoclusters = event.get(BasicCluster.class,_phnames);
+// Map<Track, SpacePoint> trposmap = (Map<Track, SpacePoint>) event.get("TrackPosMap");
+
+ int NPhoClhits = 0;
+ double EPhoClus = 0;
+
+ // check found photon clusters
+ for (BasicCluster phoclus : phoclusters)
+ {
+ int nphclhits = phoclus.getSize();
+ NPhoClhits += nphclhits;
+ EPhoClus += phoclus.getEnergy();
+ int nmchits = 0;
+ int nphomchits = 0;
+ int nneumchits = 0;
+ List<CalorimeterHit> phohits = phoclus.getCalorimeterHits();
+ for (CalorimeterHit phhit : phohits)
+ {
+ nmchits++;
+ SimCalorimeterHit hit = (SimCalorimeterHit) phhit;
+ int mcid = hit.getMCParticle(0).getPDGID();
+ double mcch = hit.getMCParticle(0).getCharge();
+// int nmchit = hit.getMCParticleCount();
+// for (int i=0; i<nmchit; ++i)
+// {
+// int mcid = hit.getPDG(i);
+ if (mcid == 22) nphomchits++;
+ if (mcch == 0.) nneumchits++;
+// }
+// nmchits += nmchit;
+// aida.cloud1D("MC Particle Count per pho hit").fill(nmchit);
+ }
+ if (nmchits>0)
+ {
+ double dnphomchits = nphomchits;
+ double dnmchits = nmchits;
+ double dnneumchits = nneumchits;
+ double hrat = dnphomchits/dnmchits;
+ double nrat = dnneumchits/nmchits;
+ aida.cloud2D("NClus hits vs NMC hits Photon").fill(nmchits,nphomchits);
+ aida.cloud1D("MC Particle Purity per Photon Clus").fill(hrat);
+ aida.cloud1D("MC Particle Purity per Neutral Clus").fill(nrat);
+// if (nmchits>nphomchits) aida.cloud1D("MC P Count gt 1 Photons").fill(hrat);
+ }
+ }
+ // Check distance between tracks at ECAL and cheated photon clusters (diag only)
+// for (Track itr : trposmap.keySet())
+// {
+// SpacePoint trsp = trposmap.get(itr);
+// double trth = trsp.theta();
+// if (trth<0) trth+= Math.PI;
+// double trph = trsp.phi();
+// if (trph<0) trph+=2*Math.PI;
+// aida.cloud2D("Theta vs Phi of Track in PhoComp").fill(trph,trth);
+ // now get cheated clusters to compare to track
+// for (CheatCluster phclus : chclusters)
+// {
+// if (phclus.getMCParticle().getPDGID() == 22)
+// {
+// double phclp[] = phclus.getPosition();
+// double phclx = phclp[0];
+// double phcly = phclp[1];
+// double phclz = phclp[2];
+// double phclR = Math.sqrt(phclx*phclx+phcly*phcly);
+// double phclth = Math.atan(phclR/phclz);
+// if (phclth<0) phclth+=Math.PI;
+// double phclph = Math.atan2(phcly,phclx);
+// if (phclph<0) phclph+=2*Math.PI;
+// aida.cloud2D("Theta vs Phi of Cheated PhoCluster").fill(phclph,phclth);
+// double distth = Math.abs(phclth-trth);
+// double distph = Math.abs(phclph-trph);
+// if (distph>Math.PI) distph=2*Math.PI-distph;
+// aida.cloud1D("Theta distance Track Cheated Pho Cluster").fill(distth);
+// aida.cloud1D("Phi distance Track Cheated Pho Cluster").fill(distph);
+// double disttrcl = Math.sqrt(distth*distth+distph*distph);
+// aida.cloud1D("Distance Track Cheated Pho Clus").fill(disttrcl);
+// }
+// }
+// }
+
+ // Check the cheated clusters also
+ double chclusE = 0;
+ int NPerfPhoCl = 0;
+ int NPhoChClhits = 0;
+ double PerfPhoClE = 0;
+ double MCPhoE = 0;
+ int NCharClus = 0;
+ for (CheatCluster chclus : chclusters)
+ {
+ chclusE += chclus.getEnergy();
+ if (phoD) aida.cloud1D("Num hits in Perfect Photon Cluster").fill(chclus.getSize());
+ if (phoD) aida.cloud2D("Num hits vs E PerfPho").fill(chclus.getEnergy(),chclus.getSize());
+ if (phoD) aida.cloud1D("Num hits over E PerfPho").fill(chclus.getSize()/chclus.getEnergy());
+ if (chclus.getMCParticle().getCharge() == 0 && chclus.getMCParticle().getPDGID() == 22)
+ {
+ NPerfPhoCl++;
+ PerfPhoClE += chclus.getEnergy();
+ MCPhoE += chclus.getMCParticle().getEnergy();
+ NPhoChClhits += chclus.getSize();
+ }
+ if (chclus.getMCParticle().getCharge() != 0)
+ {
+ NCharClus++;
+ double chclchpE = chclus.getEnergy();
+ double chclchpP = chclus.getMCParticle().getMomentum().magnitude();
+ if (phoD) aida.cloud1D("Cheated Cluster Charge E over P").fill(chclchpE/chclchpP);
+ }
+ }
+ if (chclusters.size()==1)
+ {
+ aida.cloud1D("Cheated Cluster E single chClus").fill(chclusters.get(0).getEnergy());
+ aida.cloud1D("Photon Cluster E single chClus").fill(EPhoClus);
+ if (MCPhoE>0.1)
+ {
+ aida.cloud1D("Diff Cheated Cl E MC E single chClus").fill((chclusters.get(0).getEnergy()-MCPhoE)/Math.sqrt(MCPhoE));
+ aida.cloud1D("Diff Pho ClE MC E single chClus").fill((EPhoClus-MCPhoE)/Math.sqrt(MCPhoE));
+ }
+ }
+ double dNPhoClhits = NPhoClhits;
+ double dNPhoChClhits = NPhoChClhits;
+ if (NPhoChClhits>10 && PerfPhoClE>.25)
+ {
+ double heff = dNPhoClhits/dNPhoChClhits;
+ aida.cloud1D("PhoClus Hit Efficiency").fill(heff);
+ aida.cloud1D("PhoClus E Efficiency").fill(EPhoClus/PerfPhoClE);
+// if (NPhoChClhits>NPhoClhits) aida.cloud1D("PhoClus Hit Efficiency").fill(heff);
+// if (NPhoChClhits==NPhoClhits) aida.cloud1D("PhoClus Hit Efficiency").fill(heff);
+ if (NPhoChClhits<NPhoClhits) aida.cloud1D("PhoClus Extra Hits").fill(heff);
+// if (NPhoChClhits>NPhoClhits) aida.cloud1D("PhoClus E Efficiency").fill(EPhoClus/PerfPhoClE);
+// if (NPhoChClhits==NPhoClhits) aida.cloud1D("PhoClus E Efficiency").fill(EPhoClus/PerfPhoClE);
+ aida.cloud1D("Cheated Photon Cluster E").fill(PerfPhoClE);
+ aida.cloud1D("Diff Perf Pho Clus E and MC Photon E").fill((PerfPhoClE-MCPhoE)/Math.sqrt(MCPhoE));
+ aida.cloud1D("Diff Pho Clus E and MC Pho E").fill((EPhoClus-MCPhoE)/Math.sqrt(MCPhoE));
+ aida.cloud1D("Diff Pho CLus E and Perf Pho Clus E").fill(EPhoClus-PerfPhoClE);
+ }
+ }
+
+ public void setPhoClusNames(String phnames)
+ {
+ _phnames = phnames;
+ }
+
+}