lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat
diff -N CheatReconDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CheatReconDriver.java 28 Apr 2007 00:35:09 -0000 1.1
@@ -0,0 +1,440 @@
+/*
+ * CheatReconDriver.java
+ *
+ * Created on April 27, 2007, 6:12 AM
+ *
+ * Driver to do a cheat reconstruction. Output is lists of FSparticles, clusters,
+ * tracks, and reconstructed particles
+ */
+
+package org.lcsim.contrib.Cassell.recon.Cheat;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.base.BaseTrackMC;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.recon.cluster.cheat.CheatClusterDriver;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.tracking.cheat.CheatTrackDriver;
+import org.lcsim.recon.particle.CheatParticleDriver;
+import org.lcsim.event.util.CreateFinalStateMCParticleList;
+import org.lcsim.digisim.DigiSimDriver;
+import org.lcsim.digisim.SimCalorimeterHitsDriver;
+import org.lcsim.recon.cluster.util.CalHitMapDriver;
+import java.util.HashMap;
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Map;
+import org.lcsim.util.lcio.LCIOConstants;
+import hep.physics.vec.Hep3Vector;
+/**
+ *
+ * @author cassell
+ */
+public class CheatReconDriver extends Driver
+{
+ String FSname = "GenFinalStateParticles";
+ String Tname = "RefinedCheatTracks";
+ String Cname = "RefinedCheatClusters";
+ String outRname = "CheatReconReconstructedParticles";
+ String outFSname = "CheatReconFinalStateParticles";
+ String outTname = "CheatReconTracks";
+ String outCname = "CheatReconClusters";
+ String[] diginames = {"EcalBarrDigiHits", "EcalEndcapDigiHits",
+ "HcalBarrDigiHits", "HcalEndcapDigiHits"};
+ double mintfromin = 4;
+ double mintfromout = 3;
+ double mintlength = 10.;
+ double ptcut = .1;
+ double emR = 1270.;
+ double emZ = 1680.;
+ double tccut = 200.;
+ double minneutlen = 200.;
+ Map<MCParticle,Track> tmap;
+ CheatParticleDriver cpd;
+ /** Creates a new instance of CheatSimReconParticleDriver */
+ public CheatReconDriver()
+ {
+//
+// Make a list of MCParticles to be considered
+//
+ CreateFinalStateMCParticleList cfs = new CreateFinalStateMCParticleList("Gen");
+ cfs.setCollectionName(FSname);
+ add(cfs);
+//
+// Make TrackerHit Clusters
+//
+ add(new CheatTrackDriver());
+//
+// Make CalorimeterHit Cheat Clusters
+//
+ add(new CalHitMapDriver());
+ DigiSimDriver digi = new DigiSimDriver();
+ add(digi);
+ add(new SimCalorimeterHitsDriver());
+ // output of DigiSim is hit collections - EcalBarrDigiHits, EcalEndcapDigiHits, HcalBarrDigiHits, HcalEndcapDigiHits
+ add(new CheatClusterDriver(diginames, Cname));
+ cpd = new CheatParticleDriver("dummy");
+ }
+ public void process(EventHeader event)
+ {
+ super.process(event);
+ List<MCParticle> initfs = event.get(MCParticle.class,FSname);
+ List<Track> inittracks = event.get(Track.class,Tname);
+ List<Cluster> initclusters = event.get(Cluster.class,Cname);
+ List<MCParticle> add = new ArrayList<MCParticle>();
+ List<MCParticle> remove = new ArrayList<MCParticle>();
+ List<MCParticle> tempfs = new ArrayList<MCParticle>();
+ List<MCParticle> untracked = new ArrayList<MCParticle>();
+ List<MCParticle> reconfs = new ArrayList<MCParticle>();
+ List<Track> fstracks = new ArrayList<Track>();
+ List<Cluster> fsclusters = new ArrayList<Cluster>();
+ List<Cluster> nonfsclusters = new ArrayList<Cluster>();
+ List<Cluster> untrackedchargedclusters = new ArrayList<Cluster>();
+ tmap = new HashMap<MCParticle,Track>();
+//
+// Make a map of trackable tracks
+//
+ for(Track t:inittracks)
+ {
+ MCParticle p = ((BaseTrackMC) t).getMCParticle();
+ if(isTrackable(p,t))tmap.put(p,t);
+ }
+//
+// Start with the generator final state particles
+//
+ for(MCParticle p:initfs)
+ {
+//
+// If the particle has already been removed skip it
+//
+ if(remove.contains(p))continue;
+ Hep3Vector vtx = p.getOrigin();
+ Hep3Vector ep = p.getEndPoint();
+//
+// Check if we really want the parent
+//
+ MCParticle tp = getTrackableParent(p);
+//
+// No trackable parentage. Only keep parent if really close to Cal
+//
+ if(tp == null)
+ {
+//
+// If particle is trackable, we aren't going to remove it for parentage
+//
+ if(!tmap.containsKey(p))
+ {
+ if(tooClose(vtx))
+ {
+//
+// Vertex is close to Ecal. Consider replacing with a neutral parent
+//
+ MCParticle par = getLongestNeutralParent(p);
+ if(par != null)
+ {
+ if(!hasTrackableDaughters(par))
+ {
+ add.add(par);
+ List<MCParticle> dl = new ArrayList<MCParticle>();
+ for(MCParticle d:par.getDaughters())
+ {
+ dl.add(d);
+ }
+ while(dl.size() > 0)
+ {
+ remove.add(dl.get(0));
+ for(MCParticle d:dl.get(0).getDaughters())
+ {
+ dl.add(d);
+ }
+ dl.remove(0);
+ }
+ continue;
+ }
+ }
+ }
+ }
+ }
+//
+// Have trackable parentage. Only keep parent if no trackable offspring
+// of trackable parent.
+//
+ else
+ {
+ if(!hasTrackableDaughters(tp))
+ {
+//
+// No trackable offspring, so keep trackable parent and remove all
+// offspring.
+//
+ add.add(tp);
+ List<MCParticle> dl = new ArrayList<MCParticle>();
+ for(MCParticle d:tp.getDaughters())
+ {
+ dl.add(d);
+ }
+ while(dl.size() > 0)
+ {
+ remove.add(dl.get(0));
+ for(MCParticle d:dl.get(0).getDaughters())
+ {
+ dl.add(d);
+ }
+ dl.remove(0);
+ }
+ continue;
+ }
+ }
+//
+// Done checking for parentage replacement. Now check for offspring replacement
+//
+ if(hasTrackableDaughters(p))
+ {
+ remove.add(p);
+ for(MCParticle d:p.getDaughters())
+ {
+ if(d.getSimulatorStatus().isBackscatter())continue;
+ if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
+ add.add(d);
+ }
+ }
+ else
+ {
+ if(tooClose(ep))continue;
+ if(tmap.containsKey(p))continue;
+ remove.add(p);
+ for(MCParticle d:p.getDaughters())
+ {
+ if(d.getSimulatorStatus().isBackscatter())continue;
+ if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
+ add.add(d);
+ }
+ }
+ }
+//
+// Done with pass over original FS particles. Now need to check added particles
+// for further offspring replacement
+//
+ List<MCParticle> newadd = new ArrayList<MCParticle>(add);
+ while(newadd.size() > 0)
+ {
+ List<MCParticle> add2 = new ArrayList<MCParticle>();
+ for(MCParticle p:newadd)
+ {
+ if(hasTrackableDaughters(p))
+ {
+ for(MCParticle d:p.getDaughters())
+ {
+ if(d.getSimulatorStatus().isBackscatter())continue;
+ if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
+ add2.add(d);
+ }
+ }
+ else
+ {
+ if(tooClose(p.getEndPoint()))
+ {
+ tempfs.add(p);
+ continue;
+ }
+ if(tmap.containsKey(p))
+ {
+ tempfs.add(p);
+ continue;
+ }
+ for(MCParticle d:p.getDaughters())
+ {
+ if(d.getSimulatorStatus().isBackscatter())continue;
+ if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
+ add2.add(d);
+ }
+ }
+ }
+ newadd = new ArrayList<MCParticle>(add2);
+ }
+//
+// Now have remove and add lists for a temp fs particle list. Make it.
+//
+ for(MCParticle p:initfs)
+ {
+ if(!remove.contains(p))tempfs.add(p);
+ }
+//
+// Now make the recon lists. (Removes untrackable tracks, and invisible particles
+//
+ Map<MCParticle,Cluster> fsclmap = new HashMap<MCParticle,Cluster>();
+ for(Track t:inittracks)
+ {
+ MCParticle p = ((BaseTrackMC) t).getMCParticle();
+ if(tmap.containsKey(p))
+ {
+ if(tempfs.contains(p))
+ {
+ reconfs.add(p);
+ fstracks.add(t);
+ }
+ }
+ }
+ for(Cluster c:initclusters)
+ {
+ CheatCluster cc = (CheatCluster) c;
+ MCParticle trmc = cpd.traceit(cc.getMCParticle(),tempfs);
+ if(trmc == null)
+ {
+ nonfsclusters.add(c);
+ }
+ else
+ {
+ if(fsclmap.containsKey(trmc))
+ {
+ CheatCluster bc = (CheatCluster) fsclmap.get(trmc);
+ bc.addCluster(c);
+ }
+ else
+ {
+ CheatCluster bc = new CheatCluster(trmc);
+ bc.addCluster(c);
+ fsclmap.put(trmc, bc);
+ if(!reconfs.contains(trmc))reconfs.add(trmc);
+ }
+ }
+ }
+ List<MCParticle> rm = new ArrayList<MCParticle>();
+ for(MCParticle p:reconfs)
+ {
+ if(p.getCharge() != 0)
+ {
+ if(!tmap.containsKey(p))
+ {
+ rm.add(p);
+ untrackedchargedclusters.add(fsclmap.get(p));
+ fsclmap.remove(p);
+ }
+ }
+ }
+ for(MCParticle p:rm)
+ {
+ reconfs.remove(p);
+ }
+ event.put("ReconFSParticles",reconfs,MCParticle.class,0);
+ int flag = 1<<LCIOConstants.TRBIT_HITS;
+ event.put("ReconTracks",fstracks,BaseTrackMC.class,flag);
+ flag = 1<<LCIOConstants.CLBIT_HITS;
+ event.put("ReconClusters",new ArrayList(fsclmap.values()),CheatCluster.class,flag);
+ event.put("NonFSReconClusters",nonfsclusters,CheatCluster.class,flag);
+ event.put("UntrackedChargedClusters",untrackedchargedclusters,CheatCluster.class,flag);
+ List<ReconstructedParticle> rrl = new ArrayList<ReconstructedParticle>();
+ for(MCParticle p:reconfs)
+ {
+ MCReconstructedParticle rp = new MCReconstructedParticle(p);
+ if(tmap.containsKey(p))rp.addTrack(tmap.get(p));
+ if(fsclmap.containsKey(p))rp.addCluster(fsclmap.get(p));
+ rrl.add(rp);
+ }
+ event.put("ReconPerfectReconParticles",rrl);
+ }
+ private boolean isTrackable(MCParticle p,Track t)
+ {
+ boolean tr = false;
+ Hep3Vector m = p.getMomentum();
+ if(Math.sqrt(m.x()*m.x()+m.y()*m.y()) < ptcut)return tr;
+ Hep3Vector vtx = p.getOrigin();
+ Hep3Vector ep = p.getEndPoint();
+ double d = Math.sqrt( (ep.x()-vtx.x())*(ep.x()-vtx.x()) +
+ (ep.y()-vtx.y())*(ep.y()-vtx.y()) +
+ (ep.z()-vtx.z())*(ep.z()-vtx.z()) );
+ if(d > mintlength)
+ {
+ if( (Math.sqrt(ep.x()*ep.x()+ep.y()*ep.y())>emR)||(Math.abs(ep.z())>emZ) )
+ {
+ if(t.getTrackerHits().size() >= mintfromout)tr=true;
+ }
+ else
+ {
+ if(t.getTrackerHits().size() >= mintfromin)tr=true;
+ }
+ }
+ return tr;
+ }
+ private MCParticle getTrackableParent(MCParticle p)
+ {
+ MCParticle ret = null;
+ boolean htp = false;
+ MCParticle par = p;
+ while(!htp&&(par.getParents().size() == 1))
+ {
+ par = par.getParents().get(0);
+ if(tmap.containsKey(par))htp = true;
+ }
+ if(htp)ret = par;
+ return ret;
+ }
+ private boolean hasTrackableDaughters(MCParticle p)
+ {
+ boolean htd = false;
+ List<MCParticle> check = new ArrayList<MCParticle>();
+ for(MCParticle d:p.getDaughters())
+ {
+ if(d.getSimulatorStatus().isBackscatter())continue;
+ if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
+ check.add(d);
+ }
+ while(!htd && (check.size() > 0))
+ {
+ if(tmap.containsKey(check.get(0)))htd = true;
+ else
+ {
+ List<MCParticle> dl = check.get(0).getDaughters();
+ check.remove(0);
+ for(MCParticle d:dl)
+ {
+ if(d.getSimulatorStatus().isBackscatter())continue;
+ if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
+ check.add(d);
+ }
+ }
+ }
+ return htd;
+ }
+ private boolean tooClose(Hep3Vector v)
+ {
+ boolean tc = false;
+ double vr = Math.sqrt(v.x()*v.x() + v.y()*v.y());
+ if( (vr > emR-tccut)||(Math.abs(v.z()) > emZ-tccut) )
+ {
+ tc = true;
+ }
+ return tc;
+ }
+ private MCParticle getLongestNeutralParent(MCParticle p)
+ {
+ MCParticle ret = null;
+ double length = minneutlen;
+ MCParticle par = p;
+ while(par.getParents().size() == 1)
+ {
+ par = par.getParents().get(0);
+ if(par.getCharge() != 0)
+ {
+ Hep3Vector v = par.getOrigin();
+ Hep3Vector e = par.getEndPoint();
+ double l = Math.sqrt( (e.x()-v.x())*(e.x()-v.x()) +
+ (e.y()-v.y())*(e.y()-v.y()) +
+ (e.z()-v.z())*(e.z()-v.z()) );
+ if(l > length)
+ {
+ length = l;
+ ret = par;
+ }
+ }
+ }
+ return ret;
+ }
+
+}