lcsim/src/org/lcsim/contrib/SteveMagill

diff -N TrMipCompDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ TrMipCompDriver.java	29 May 2008 21:32:47 -0000	1.1
@@ -0,0 +1,133 @@

+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.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 TrMipCompDriver extends Driver
+{
+  public TrMipCompDriver() {
+  }
+
+  private String _tmclusname;
+  private AIDA aida = AIDA.defaultInstance();
+  private boolean perfPFAD = true;
+  
+  public void process(EventHeader event) 
+  {
+      List<CheatCluster> chclusters = event.get(CheatCluster.class,"PerfectCheatClusters");
+
+      int NMipCl = 0;
+      int nhitstm = 0;  // number of hits in all mip clusters
+      int nchhitstm = 0;  // number of charged particle hits in all mip clusters      
+      try
+      {
+      List<BasicCluster> tmclusters = event.get(BasicCluster.class,_tmclusname);
+
+      for (BasicCluster tmclus : tmclusters)
+      {
+          NMipCl++;
+          int numhits = tmclus.getSize();
+          nhitstm += tmclus.getSize();  // sum hits from all mips in event
+          double tmclE = tmclus.getEnergy();
+          aida.cloud1D("Tr Mip Clus E").fill(tmclE);
+          aida.cloud1D("Tr Mip Clus Size").fill(numhits);
+          aida.cloud2D("Tr Mip NumHits vs E").fill(tmclE,numhits);
+          List<CalorimeterHit> tmcalhits = tmclus.getCalorimeterHits();
+          int numchits = 0;
+          for (CalorimeterHit tmhit : tmcalhits)
+          {
+              SimCalorimeterHit hit = (SimCalorimeterHit) tmhit;
+              double hitch = Math.abs(hit.getMCParticle(0).getCharge());
+              if (hitch>0) numchits++;
+              if (hitch>0) nchhitstm++;
+          }
+          if (numhits>0)
+          {
+              double dnumhits = numhits;
+              double dnumchits = numchits;
+              double hrat = dnumchits/dnumhits;
+              aida.cloud2D("NClus hits vs NMC hits Tr Mips").fill(numchits,numhits);
+              aida.cloud1D("MC Particle Purity per Tr Mip Clus").fill(hrat);
+          }
+      }
+      double dnhitstm = nhitstm;
+      double dnchhitstm = nchhitstm;
+      double mipTot = dnchhitstm/dnhitstm;
+      if (nhitstm>0) aida.cloud1D("MC Particle Purity for all Tr Mip Clus").fill(mipTot);
+      
+      }     
+      catch(java.lang.IllegalArgumentException ex)
+      {
+          System.out.println("requested object not found in event " + _tmclusname);      
+      }
+      
+      if (NMipCl == 1)
+      {
+          aida.cloud1D("Num Hits Mip Clus single cluster").fill(nhitstm);
+      }
+      
+      //  Check the cheated clusters also - only good for muons 
+      double chclusE = 0;
+      int NCheatCl = 0;
+      double PerfChHClE = 0;
+      double MCChHE = 0;
+      int nhitsingcl = 0;
+      for (CheatCluster chclus : chclusters)
+      {
+          chclusE += chclus.getEnergy();
+          aida.cloud1D("Num hits in Perfect Cheat Cluster Mips").fill(chclus.getSize());
+          aida.cloud2D("Num hits vs E PerfCheatCl Mips").fill(chclus.getEnergy(),chclus.getSize());
+          aida.cloud1D("Num hits over E PerfCheatCl Mips").fill(chclus.getSize()/chclus.getEnergy());
+          if (chclus.getMCParticle().getCharge() != 0)
+          {
+              NCheatCl++;
+              nhitsingcl += chclus.getSize();
+              PerfChHClE += chclus.getEnergy();
+              MCChHE += chclus.getMCParticle().getEnergy();
+              double chclchpE = chclus.getEnergy();
+              double chclchpP = chclus.getMCParticle().getMomentum().magnitude();
+              aida.cloud1D("Cheated Cluster Charge E over P Mips").fill(chclchpE/chclchpP);
+          }
+      }
+      if (NCheatCl == 1)
+      {
+          aida.cloud1D("Cheated Cluster E single cluster Mips").fill(chclusE);
+          aida.cloud1D("Num Hits in Cheated Cluster single cluster Mips").fill(nhitsingcl);
+          if (NMipCl == 1) aida.cloud1D("Diff NHits TM Cheat Clusters").fill(nhitstm-nhitsingcl);
+      }
+  }
+  
+  public void setTrMipClusName(String tmclusname)
+  {
+      _tmclusname = tmclusname;
+  }
+  
+}