lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon/Cheat
diff -N ClusterLevelCheater.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClusterLevelCheater.java 26 Oct 2010 20:20:14 -0000 1.1
@@ -0,0 +1,220 @@
+package org.lcsim.contrib.Cassell.recon.Cheat;
+import org.lcsim.util.aida.AIDA;
+import java.util.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.contrib.Cassell.recon.analysis.*;
+import hep.physics.vec.*;
+import org.lcsim.recon.cluster.nn.NearestNeighborClusterer;
+import org.lcsim.recon.cluster.util.*;
+
+public class ClusterLevelCheater extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ int ievt;
+ String epre;
+ String hpre;
+ String[] EcalCl;
+ String[] HcalCl;
+ String inRecon;
+ String outRecon;
+ String fsp;
+ String fstr;
+ ClusterEnergyCalculator scalc;
+ ClusterEnergyCalculator calc;
+ int minCalhits = 5;
+ Clusterer clusterer;
+ Hep3Vector mom;
+ int ncuts = 2;
+ int[] cut = {2,5};
+ Map<MCParticle,List<Cluster>> emap;
+ Map<MCParticle,List<Cluster>> hmap;
+ public ClusterLevelCheater(String Eprefix,String Hprefix,String[] ecalcl,String[] hcalcl, String inrecon, String outrecon,String fsmcp,String fsmct)
+ {
+ this(Eprefix,Hprefix,ecalcl,hcalcl,inrecon,outrecon,fsmcp,fsmct,false);
+ }
+ public ClusterLevelCheater(String Eprefix,String Hprefix,String[] ecalcl,String[] hcalcl, String inrecon, String outrecon,String fsmcp,String fsmct,boolean analog)
+ {
+ ievt = 0;
+ epre = Eprefix;
+ hpre = Hprefix;
+ EcalCl = ecalcl;
+ HcalCl = hcalcl;
+ inRecon = inrecon;
+ outRecon = outrecon;
+ fsp = fsmcp;
+ fstr = fsmct;
+ scalc = new QPhotonClusterEnergyCalculator();
+ calc = new QNeutralHadronClusterEnergyCalculator(analog);
+ clusterer = new NearestNeighborClusterer(4,4,2,0,0.);
+ }
+ Map<MCParticle,List<Cluster>> getEcalMap(){return emap;}
+ Map<MCParticle,List<Cluster>> getHcalMap(){return hmap;}
+
+ protected void process(EventHeader event)
+ {
+// System.out.println("Entering ClusterLevelCheater.process");
+//
+// Get all the necessary collections from event
+//
+ List<MCParticle> fst = event.get(MCParticle.class,fstr);
+ List<MCParticle> fs = event.get(MCParticle.class,fsp);
+ List<Cluster> ecl = new ArrayList<Cluster>();
+ List<Cluster> hcl = new ArrayList<Cluster>();
+ for(int i=0;i<EcalCl.length;i++)
+ {
+ ecl.addAll(event.get(Cluster.class,EcalCl[i]));
+ }
+ for(int i=0;i<HcalCl.length;i++)
+ {
+ hcl.addAll(event.get(Cluster.class,HcalCl[i]));
+ }
+ List<ReconstructedParticle> inrl = event.get(ReconstructedParticle.class,inRecon);
+ List<ReconstructedParticle>[] outrl = new List[ncuts];
+ for(int i=0;i<ncuts;i++)
+ {
+ outrl[i] = new ArrayList<ReconstructedParticle>();
+ }
+// System.out.println("ClusterLevelCheater: making ClusterPurutyPlots");
+ ClusterPurityPlots cpp = new ClusterPurityPlots(fs,fst);
+ cpp.makePlots(epre,ecl);
+ emap = cpp.getMap();
+ cpp.makePlots(hpre,hcl);
+ hmap = cpp.getMap();
+ for(int i=0;i<ncuts;i++)
+ {
+ for(ReconstructedParticle inrp:inrl)
+ {
+ BaseReconstructedParticle orp;
+ MCParticle p = null;
+ if(inrp.getParticles().size() > 0)p = ( (MCReconstructedParticle) inrp.getParticles().get(0) ).getMCParticle();
+ if(inrp.getCharge() != 0)
+ {
+ orp = new BaseReconstructedParticle(inrp.getEnergy(),inrp.getMomentum());
+ orp.setMass(inrp.getMass());
+ orp.setCharge(inrp.getCharge());
+ orp.setReferencePoint(inrp.getReferencePoint());
+ orp.setParticleIdUsed(inrp.getParticleIDUsed());
+ orp.addParticleID(inrp.getParticleIDUsed());
+ for(Track t:inrp.getTracks()){orp.addTrack(t);}
+ int necl = 0;
+ double maxe = 0.;
+ if(p != null)
+ {
+ if(emap.containsKey(p))
+ {
+ for(Cluster c:emap.get(p))
+ {
+ if(c.getCalorimeterHits().size() > cut[i])
+ {
+ orp.addCluster(c);
+ necl++;
+ if(c.getEnergy() > maxe)maxe = c.getEnergy();
+ }
+ }
+ }
+ aida.profile1D("track/mean # Ecal clusters > "+cut[i]+" hits vs E",100,0.,100.).fill(p.getEnergy(),necl);
+ if(i == 0)aida.profile1D("track/mean max Ecal clusterE vs E",100,0.,100.).fill(p.getEnergy(),maxe);
+ int nhcl = 0;
+ maxe = 0.;
+ if(hmap.containsKey(p))
+ {
+ for(Cluster c:hmap.get(p))
+ {
+ if(c.getCalorimeterHits().size() > i)
+ {
+ orp.addCluster(c);
+ nhcl++;
+ if(c.getEnergy() > maxe)maxe = c.getEnergy();
+ }
+ }
+ }
+ aida.profile1D("track/mean # Hcal clusters > "+cut[i]+" hits vs E",100,0.,100.).fill(p.getEnergy(),nhcl);
+ if(i == 0)aida.profile1D("track/mean max Hcal clusterE vs E",100,0.,100.).fill(p.getEnergy(),maxe);
+ orp.addParticle(inrp.getParticles().get(0));
+ }
+ outrl[i].add(orp);
+ }
+ else
+ {
+ BasicCluster nc = new BasicCluster();
+ double E = 0.;
+ String ntype = "photon/";
+ if(inrp.getMass() != 0.)ntype = "n hadron/";
+ int necl = 0;
+ double maxe = 0.;
+ if(emap.containsKey(p))
+ {
+ for(Cluster c:emap.get(p))
+ {
+ if(c.getCalorimeterHits().size() > cut[i])
+ {
+ nc.addCluster(c);
+ necl++;
+ if(c.getEnergy() > maxe)maxe = c.getEnergy();
+ }
+ }
+ }
+ aida.profile1D(ntype+"mean # Ecal clusters > "+cut[i]+" hits vs E",100,0.,100.).fill(p.getEnergy(),necl);
+ if(i == 0)aida.profile1D(ntype+"mean max Ecal clusterE vs E",100,0.,100.).fill(p.getEnergy(),maxe);
+ int nhcl = 0;
+ maxe = 0.;
+ if(hmap.containsKey(p))
+ {
+ for(Cluster c:hmap.get(p))
+ {
+ if(c.getCalorimeterHits().size() > cut[i])
+ {
+ nc.addCluster(c);
+ nhcl++;
+ if(c.getEnergy() > maxe)maxe = c.getEnergy();
+ }
+ }
+ }
+ aida.profile1D(ntype+"mean # Hcal clusters > "+cut[i]+" hits vs E",100,0.,100.).fill(p.getEnergy(),nhcl);
+ if(i == 0)aida.profile1D(ntype+"mean max Hcal clusterE vs E",100,0.,100.).fill(p.getEnergy(),maxe);
+ if(nc.getCalorimeterHits().size() >= minCalhits)
+ {
+ if(inrp.getMass() == 0.)E = scalc.getEnergy(nc);
+ else E = calc.getEnergy(nc);
+ if(E > inrp.getMass())
+ {
+ List<Cluster> clist = clusterer.createClusters(nc.getCalorimeterHits());
+ double[] pos = {0.,0.,0.};
+ int maxhits = 0;
+ for(Cluster ccc:clist)
+ {
+ if(ccc.getCalorimeterHits().size() > maxhits)
+ {
+ pos = ccc.getPosition();
+ maxhits = ccc.getCalorimeterHits().size();
+ }
+ }
+ double norm = Math.sqrt(E*E - inrp.getMass()*inrp.getMass())/
+ Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
+ mom = new BasicHep3Vector(pos[0]*norm, pos[1]*norm, pos[2]*norm);
+ orp = new BaseReconstructedParticle(E,mom);
+ orp.setMass(inrp.getMass());
+ orp.setCharge(inrp.getCharge());
+ orp.setReferencePoint(inrp.getReferencePoint());
+ orp.setParticleIdUsed(inrp.getParticleIDUsed());
+ orp.addParticleID(inrp.getParticleIDUsed());
+ for(Track t:inrp.getTracks()){orp.addTrack(t);}
+ for(Cluster c:nc.getClusters()){orp.addCluster(c);}
+ outrl[i].add(orp);
+ orp.addParticle(inrp.getParticles().get(0));
+ }
+ }
+ }
+ }
+ event.put(outRecon+"cut"+cut[i],outrl[i]);
+ }
+ ievt++;
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N DebugReconDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DebugReconDriver.java 26 Oct 2010 20:20:14 -0000 1.1
@@ -0,0 +1,82 @@
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.contrib.Cassell.recon.Cheat.ClusterLevelCheater;
+import org.lcsim.recon.cheater.*;
+import org.lcsim.util.Driver;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.pfa.structural.SetUpPFA;
+import org.lcsim.recon.ui.ReconDriver;
+
+/**
+ *
+ * @author cassell
+ */
+public class DebugReconDriver extends Driver
+{
+ String CheatReconRname = "ReconPerfectReconParticles";
+ String PPRPflowRname = "PPRReconParticles";
+ String CheatReconFSname = "ReconFSParticles";
+ String ExtraTracks = "ExtraTracks";
+ String CheatReconFSTrackedname = "TrackedReconFSParticles";
+ String CheatCname = "ReconClusters";
+ String[] EcalNames = {"EcalBarrDigiHits", "EcalEndcapDigiHits"};
+ String[] HcalNames = {"HcalBarrDigiHits", "HcalEndcapDigiHits"};
+ String[] hitNames = {EcalNames[0],EcalNames[1],HcalNames[0],HcalNames[1]};
+ String clusterername = "DTreeClusters";
+ String[] OEClNames = {"EcalBarrDigiHits"+clusterername, "EcalEndcapDigiHits"+clusterername};
+ String[] OHClNames = {"HcalBarrDigiHits"+clusterername, "HcalEndcapDigiHits"+clusterername};
+ String[] REClNames = {"EcalBarrDigiHitsReclustered"+clusterername, "EcalEndcapDigiHitsReclustered"+clusterername};
+ String[] RHClNames = {"HcalBarrDigiHitsReclustered"+clusterername, "HcalEndcapDigiHitsReclustered"+clusterername};
+ String ClLevelPFlowRname = clusterername+"ReconParticles";
+ String[] prefix = {"PPR",clusterername+">2hits",
+ clusterername+">5hits"};
+ String[] rnames = {PPRPflowRname,ClLevelPFlowRname+"cut2",
+ ClLevelPFlowRname+"cut5"};
+ private boolean useNewInitialMipFinding = false;
+ SetUpPFA setup;
+ public void setUseNewInitialMipFinding(boolean x)
+ {
+ useNewInitialMipFinding = x;
+ setup.setUseNewInitialMipFinding(useNewInitialMipFinding);
+ }
+ public DebugReconDriver()
+ {
+ ReconDriver rd = new ReconDriver();
+ rd.setUseNewInitialMipFinding(true);
+ add(rd);
+//
+// Do the cheating reconstruction (includes Digisim)
+//
+ CheatReconDriverRealTracking crd = new CheatReconDriverRealTracking();
+ crd.setCheatReconstructedParticleOutputName(CheatReconRname);
+ crd.setCheatFSParticleOutputName(CheatReconFSname);
+ add(crd);
+
+
+
+//
+// Make the perfect pattern recognition pflow reconstructed particles from the
+// Cheat Recon particles
+//
+ PPRDriverRealTracking pprd = new PPRDriverRealTracking(CheatReconRname,PPRPflowRname,ExtraTracks);
+ ClusterEnergyCalculator pcec = new QPhotonClusterEnergyCalculator();
+ pprd.setPhotonClusterEnergyCalculator(pcec);
+ ClusterEnergyCalculator nhcec = new QNeutralHadronClusterEnergyCalculator();
+ pprd.setNeutralHadronClusterEnergyCalculator(nhcec);
+ add(pprd);
+//
+// Recluster the Ecal
+//
+ add(new CoreReclusterDriver(OEClNames[0],REClNames[0]));
+ add(new CoreReclusterDriver(OEClNames[1],REClNames[1]));
+//
+// Recluster the Hcal
+//
+ add(new DigitalCoreReclustererDriver(OHClNames[0],RHClNames[0]));
+ add(new DigitalCoreReclustererDriver(OHClNames[1],RHClNames[1]));
+//
+// Add cluster level cheater
+//
+ add(new ClusterLevelCheater("E"+clusterername+"recl","H"+clusterername,REClNames,RHClNames,PPRPflowRname,ClLevelPFlowRname,
+ CheatReconFSname,CheatReconFSTrackedname));
+ }
+}