lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -N ConnectTheDotsConverter.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ConnectTheDotsConverter.java 8 Dec 2007 00:23:31 -0000 1.1
@@ -0,0 +1,484 @@
+/*
+ * ConnectTheDotsConverter.java
+ *
+ * Created on October 18, 2007, 11:53 PM
+ *
+ * To change this template, choose Tools | Template Manager
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.CosminDeaconu;
+
+import hep.graphics.heprep.HepRepFactory;
+import hep.graphics.heprep.HepRepInstance;
+import hep.graphics.heprep.HepRepInstanceTree;
+import hep.graphics.heprep.HepRepType;
+import hep.graphics.heprep.HepRepTypeTree;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.awt.Color;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.List;
+import org.lcsim.contrib.tracking.TrackerHitCheater;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.util.heprep.HepRepCollectionConverter;
+import org.lcsim.util.heprep.LCSimHepRepConverter;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.event.MCParticle.SimulatorStatus;
+
+/**
+ * DESCRIPTION:
+ *
+ * This converter is meant to more realistically draw MCParticle trajectories in
+ * the event display. Rather than relying on initial data for the MCParticle as
+ * the default MCParticleConverter does, this converter "connects the dots"
+ * between SimTrackerHits belonging to the same MCParticle.
+ *
+ *
+ * USAGE:
+ *
+ * This converter will operate on SimTrackerHitCollections called "AllSimHits"
+ * (or whatever you change the variable COLLECTION_NAME to). You can use the
+ * PutAllHitsInOneCollection class in my contrib folder to make this collection
+ * for you.
+
+ * If for some reason you want lines in between the SimTrackerHits instead of
+ * helices, set FORCE_LINE_DRAW to true.
+ *
+ * You may increase the resolution of the helices by modifying D_ALPHA...
+ *
+ *
+ * "INSTALLATION":
+ *
+ * To enable this converter, you must modify
+ * org.lcsim.util.heprep.LCSimHepRepConverter.java in the following way:
+ *
+ * 1) Add "org.lcsim.contrib.CosminDeaconu.ConnectTheDotsConverter"
+ * to the imports
+ *
+ * 2) Add "register(new ConnectTheDotsConverter());" to the first try
+ * statement.
+ *
+ * It should look something like this:
+ *
+ * try
+ * {
+ * factory = HepRepFactory.create();
+ * register(new CalorimeterHitConverter());
+ * register(new ClusterConverter());
+ *
+ * ...
+ *
+ * register(new ConnectTheDotsConverter());
+ * }
+ *
+ * Note that the converter might fail if SimTrackerHits don't have momenta
+ * associated with them. I believe that this feature is enabled by default,
+ * but if you've changed THBIT_MOMENTUM in LCIOConstants, you might have
+ * problems. In that case, you'll probably have to set FORCE_LINE_DRAW to true
+ * and you'll only get lines connecting the SimTrackerHits...
+ *
+ * ISSUES:
+ *
+ * 1) Helices may appear "cut-off" near a hit, increasing resolution may or
+ * may not help here.
+ *
+ * 2) Be aware that anything with over 3000 points is split into multiple. This
+ * is due to a hard-coded limitation in Wired that limits the nubmer of points
+ * an instance may have. If anyone has all the necessary dependencies to compile
+ * the Wired4 plugin, please send them to me =D.
+ *
+ *
+ * @author Cosmin Deaconu
+ *
+ */
+public class ConnectTheDotsConverter implements HepRepCollectionConverter{
+
+
+ private static final double D_ALPHA=5; //5 mm. Make sure to change MAX_SEX_PER_D appropriately too if you change this
+ private static final double MIN_ENERGY = 0; //allows cut on energy
+ private static final String COLLECTION_NAME = "AllSimHits";
+ private static final String IDENTIFIER_STRING = "ConnectTheDotsMCParticles";
+ private static final int MAX_SEGS_PER_D = 10;//maximum segments per Euclidean distance before just drawing a line...
+ // this is sort of an upper limit on complexity and also a safeguard
+ // against tracks with infinite points...
+ private static final boolean FORCE_STOP_AT_TRACKING_REGION_BOUNDARY = false;
+
+ private static final Color ELECTRON_COLOR = Color.PINK; //color for e+/e-
+ private static final Color PION_COLOR = Color.RED; //color for pi+/pi-
+ private static final Color MUON_COLOR = Color.CYAN; //color for mu+/mu-
+ private static final Color PROTON_COLOR = Color.LIGHT_GRAY; //color for p, p_bar
+ private static final Color KAON_COLOR = Color.YELLOW; //color for K+,K-
+ private static final Color OTHER_COLOR = Color.MAGENTA; //color for anything else
+
+ private static final int MAX_POINTS = 3000; //this is the max number of points per instance in Wired4
+
+ //always draw straight connecting lines, in case you don't have momentum info
+ //in SimTrackerHits... or if you really like straight lines
+ private static final boolean FORCE_LINE_DRAW = false;
+
+ private HepRepType type;
+ private HepRepType electronType;
+ private HepRepType pionType;
+ private HepRepType muonType;
+ private HepRepType protonType;
+ private HepRepType kaonType;
+ private HepRepType genericType;
+
+ /** Creates a new instance of ConnectTheDotsConverter */
+ public boolean canHandle(Class k) {
+ return SimTrackerHit.class.isAssignableFrom(k);
+ }
+
+ public void convert(EventHeader event, List collection, HepRepFactory factory, HepRepTypeTree typeTree, HepRepInstanceTree instanceTree) {
+
+ LCMetaData meta = event.getMetaData(collection);
+
+ if (meta.getName()!=COLLECTION_NAME) return; //only process events with the right name
+
+
+ // set up the instance types
+ type = factory.createHepRepType(typeTree, IDENTIFIER_STRING);
+
+ type.addAttValue("layer",LCSimHepRepConverter.PARTICLES_LAYER);
+ type.addAttValue("drawAs","Line");
+ type.addAttDef("type","Particle Type", "physics", "");
+ type.addAttDef("energy","Initial Particle Energy","physics","GeV");
+ type.addAttDef("p","Initial Particle Energy","physics","GeV");
+ type.addAttDef("numhits","Number of Hits","physics","");
+ type.addAttDef("pt","Initial Particle Transverse Energy","physics","GeV");
+ type.addAttDef("time","Particle Production Time","physics","nanoseconds");
+ type.addAttDef("isCurler","Whether or not particle is a curler, by some definition","physics","");
+ type.addAttDef("particleStatus","the Simulator Status of the particle","physics","");
+ type.addAttDef("particleOriginRadius","The radius in mm from the ip to the origin of the MCParticle","physics","");
+ type.addAttDef("particleOriginZ","Z coordinate of particle origin","physics","");
+ type.addAttDef("split","path split because of too many points","!physics","");
+
+ electronType = factory.createHepRepType(type,"Electron");
+ electronType.addAttValue("color",ELECTRON_COLOR);
+
+ pionType = factory.createHepRepType(type,"Pion");
+ pionType.addAttValue("color",PION_COLOR);
+
+ muonType = factory.createHepRepType(type,"Muon");
+ muonType.addAttValue("color",MUON_COLOR);
+
+ protonType = factory.createHepRepType(type,"Proton");
+ protonType.addAttValue("color",PROTON_COLOR);
+
+ kaonType = factory.createHepRepType(type,"Kaon");
+ kaonType.addAttValue("color",KAON_COLOR);
+
+ genericType = factory.createHepRepType(type,"Other Particle");
+ genericType.addAttValue("color",OTHER_COLOR);
+
+ //We need to know the particles belonging to each hit...
+ HashMap<MCParticle,ArrayList<SimTrackerHit>> hits_by_particle = new HashMap<MCParticle,ArrayList<SimTrackerHit>>();
+
+ for (SimTrackerHit hit : (List<SimTrackerHit>) collection) {
+ MCParticle p = hit.getMCParticle();
+ if (hits_by_particle.containsKey(p)) {
+ ArrayList<SimTrackerHit> hitlist =hits_by_particle.get(p);
+ if (!hitlist.contains(hit)) hitlist.add(hit);
+ hits_by_particle.put(p,hitlist);
+ }
+
+ else {
+ ArrayList<SimTrackerHit> hitlist = new ArrayList<SimTrackerHit>();
+ hitlist.add(hit);
+ hits_by_particle.put(p,hitlist);
+ }
+ }
+
+ //get Helix for swimming
+ HelixSwimmer helix = new HelixSwimmer(event.getDetector().getFieldMap().getField(new double[]{0,0,0})[2]);
+
+ // loop through each particle
+ for (MCParticle p : hits_by_particle.keySet()) {
+ if (p.getEnergy()<MIN_ENERGY) continue;
+
+ ArrayList<SimTrackerHit> hitlist = hits_by_particle.get(p); //get hits from this particle
+
+ TimeCompare<SimTrackerHit> cmp = new TimeCompare<SimTrackerHit>();
+ Collections.sort(hitlist,cmp); //sort hits by time
+
+ InstanceContext context = new InstanceContext(p,hitlist,factory,instanceTree); //convenience class
+ HepRepInstance instance = getInstance(context); //create the instance
+
+ //loop through simhits
+ for (int i = 0; i< hitlist.size(); i++) {
+ SimTrackerHit h = hitlist.get(i);
+
+ //if flag set, connect the hits with lines
+ if (FORCE_LINE_DRAW) {
+
+ //draw the first point if first iteration
+ if (i==0) instance = makePoint(context,instance,p.getOrigin().x(),p.getOrigin().y(),p.getOrigin().z());
+ instance = makePoint(context,instance,h.getPoint()[0],h.getPoint()[1],h.getPoint()[2]);
+
+ //if last point, try to draw line to final point
+ if (i==hitlist.size()-1) {
+ try {
+ double[] end = p.getEndPoint().v();
+ instance = makePoint(context,instance,end[0],end[1],end[2]);
+ }
+ catch(RuntimeException re) {
+ //for now, this does nothing... presumably we'd want to draw a line to the tracker boundary or something
+ }
+ }
+ }
+
+ //otherwise, connect hits with helices
+ else {
+
+ double zStop; //swim the helix till we get here in z...
+ double charge = p.getCharge();
+
+ //case where we're not dealing with the last hit
+ if (i<hitlist.size()-1) {
+
+ if (i==0) {//draw from origin to first point if it's the first point
+ helix.setTrack(p.getMomentum(),new SpacePoint(p.getOrigin()),(int) charge);
+ zStop = hitlist.get(0).getPoint()[2];
+
+ int segs =1+(int) (helix.getDistanceToZ(zStop)/D_ALPHA);
+ double lDist = VecOp.sub(new BasicHep3Vector(hitlist.get(0).getPoint()),p.getOrigin()).magnitude();
+
+ //if segs exceeds MAX_SEGS...we just draw a straight line...
+ if (segs<MAX_SEGS_PER_D*lDist && segs>0) {
+ for (int j = 0; j<=segs; j++) {
+ SpacePoint pt = helix.getPointAtLength(D_ALPHA*j);
+ instance = makePoint(context,instance,pt.x(),pt.y(),pt.z());
+ }
+ }
+ else {
+ instance = makePoint(context,instance,p.getOriginX(),p.getOriginY(),p.getOriginZ());
+ instance = makePoint(context,instance,h.getPoint()[0],h.getPoint()[1],h.getPoint()[2]);
+ }
+ }
+
+ zStop = hitlist.get(i+1).getPoint()[2];
+
+ Hep3Vector momentum = new BasicHep3Vector(h.getMomentum());
+ SpacePoint start = new SpacePoint(new BasicHep3Vector(h.getPoint()));
+
+ helix.setTrack(momentum,start,(int)charge);
+
+ int segs = 1+(int) (helix.getDistanceToZ(zStop)/D_ALPHA);
+ double lDist = VecOp.sub(new BasicHep3Vector(hitlist.get(i+1).getPoint()),start).magnitude();
+
+ if (segs<MAX_SEGS_PER_D*lDist) {
+ for (int j = 0; j<=segs; j++) {
+ SpacePoint pt = helix.getPointAtLength(D_ALPHA*j);
+ instance = makePoint(context,instance,pt.x(),pt.y(),pt.z());
+ }
+ }
+ else {
+ instance = makePoint(context,instance,h.getPoint()[0],h.getPoint()[1],h.getPoint()[2]);
+ instance = makePoint(context,instance,hitlist.get(i+1).getPoint()[0],hitlist.get(i+1).getPoint()[1],hitlist.get(i+1).getPoint()[2]);
+ }
+ }
+
+ //case where we ARE dealing with the last hit
+ else {
+ Hep3Vector momentum = new BasicHep3Vector(h.getMomentum());
+ SpacePoint start = new SpacePoint(new BasicHep3Vector(h.getPoint()));
+ helix.setTrack(momentum,start,(int)charge);
+
+ double r = event.getDetector().getConstants().get("tracking_region_radius").getValue();
+ double z = event.getDetector().getConstants().get("tracking_region_zmax").getValue();
+
+ try{
+ zStop = p.getEndPoint().z(); //this might thrown an exception if there is no endpoint defined
+
+ int segs = (int) Math.abs(helix.getDistanceToZ(zStop)/D_ALPHA);
+ double lDist = VecOp.sub(start,p.getEndPoint()).magnitude();
+
+ if (segs<MAX_SEGS_PER_D*lDist && segs>0) {
+ for (int j = 0; j<=segs; j++) {
+ SpacePoint pt = helix.getPointAtLength(D_ALPHA*j);
+ instance = makePoint(context,instance,pt.x(),pt.y(),pt.z());
+
+ if (FORCE_STOP_AT_TRACKING_REGION_BOUNDARY) {if (pt.rxy()>r || Math.abs(pt.z()) > z) break;}
+ }
+ }
+
+ else {
+ instance = makePoint(context,instance,h.getPoint()[0],h.getPoint()[1],h.getPoint()[2]);
+ instance = makePoint(context,instance,p.getEndPoint().x(),p.getEndPoint().y(),p.getEndPoint().z());
+ }
+ }
+
+ //if no endpoint...
+ catch (RuntimeException re) {
+ int j = 0;
+ int segs = 0;
+ while (true) {
+ SpacePoint pt = helix.getPointAtLength(D_ALPHA*j);
+ double lDist = VecOp.sub(new BasicHep3Vector(h.getPoint()),new BasicHep3Vector(pt.v())).magnitudeSquared();
+
+ if (segs>MAX_SEGS_PER_D*lDist) break;
+ if (pt.rxy()>r || Math.abs(pt.z()) > z) break;
+
+ instance = makePoint(context,instance,pt.x(),pt.y(),pt.z());
+ j++;
+ segs ++;
+ }
+ }
+ }//else last hit
+ }//else draws helices
+ }// loop through simhits
+ }//loop through particles
+ }//end convert()
+
+
+ //int so you can cut on it...
+ int isCurler(ArrayList<SimTrackerHit> hitlist) {
+
+ TrackerHitCheater thcheater = new TrackerHitCheater();
+ int numhits = thcheater.makeTrackerHits(hitlist).size();
+ ArrayList<Layer> layershit = new ArrayList<Layer>();
+
+ for (SimTrackerHit h : hitlist) {
+ Layer l = new Layer(h.getLayer(),h.getSubdetector().getName());
+ if (!layershit.contains(l)) layershit.add(l);
+ }
+
+ if (numhits>layershit.size()) return 1;
+ return 0;
+ }
+
+ //borrowed/modified from MCParticleTableModel
+ String getSimulatorString(SimulatorStatus status) {
+ List<String> s = new ArrayList<String>();
+ if (status.hasLeftDetector()) s.add("Left");
+ if (status.isBackscatter()) s.add("Backscatter");
+ if (status.isCreatedInSimulation()) s.add("Created In Simulation");
+ if (status.isDecayedInCalorimeter()) s.add("Decayed In Calorimeter");
+ if (status.isDecayedInTracker()) s.add("Decayed In Tracker");
+ if (status.vertexIsNotEndpointOfParent()) s.add("Vertex not endpoint of parent");
+ if (status.isStopped()) s.add("Stopped");
+ StringBuffer buff = new StringBuffer();
+ for (int i=0; i<s.size(); i++) buff.append(s.get(i)+",");
+ buff.setLength(Math.max(0,buff.length()-1));
+ return buff.toString();
+ }
+
+ //wrapper class around createHepRepPoint that will create a new instance if too many points. Returns the instance for the next points to build on.
+ private HepRepInstance makePoint(InstanceContext context, HepRepInstance instance, double x, double y, double z) {
+
+ if (instance.getPoints().size()<MAX_POINTS) {
+ context.factory.createHepRepPoint(instance,x,y,z);
+ return instance;
+ }
+
+ else {
+ HepRepInstance instance2 = getInstance(context);
+ context.factory.createHepRepPoint(instance2,x,y,z);
+ instance.addAttValue("split",1);
+ instance2.addAttValue("split",1);
+ return instance2;
+ }
+ }
+
+
+ //creates a new instance from a particle and its hitlist
+ private HepRepInstance getInstance(InstanceContext c) {
+ String ptype = c.p.getType().getName();
+ HepRepInstance instance;
+
+ if (ptype.equals("e-")||ptype.equals("e+")) instance = c.factory.createHepRepInstance(c.instanceTree,electronType);
+ else if (ptype.equals("mu-")||ptype.equals("mu+")) instance = c.factory.createHepRepInstance(c.instanceTree,muonType);
+ else if (ptype.equals("pi-")||ptype.equals("pi+")) instance = c.factory.createHepRepInstance(c.instanceTree,pionType);
+ else if (ptype.equals("p")||ptype.equals("p_bar")) instance = c.factory.createHepRepInstance(c.instanceTree,protonType);
+ else if (ptype.equals("K+")||ptype.equals("K-")) instance = c.factory.createHepRepInstance(c.instanceTree,kaonType);
+ else instance = c.factory.createHepRepInstance(c.instanceTree,genericType);
+
+ double x = c.p.getMomentum().x();
+ double y = c.p.getMomentum().y();
+ double pT = Math.sqrt(x*x + y*y);
+
+ instance.addAttValue("type",ptype);
+ instance.addAttValue("energy",c.p.getEnergy());
+ instance.addAttValue("pt",pT);
+ instance.addAttValue("p",c.p.getMomentum().magnitude());
+ instance.addAttValue("numhits",c.hitlist.size());
+ instance.addAttValue("time",c.p.getProductionTime());
+ instance.addAttValue("hashcode",c.p.hashCode());
+ instance.addAttValue("isCurler",isCurler(c.hitlist));
+ instance.addAttValue("particleOriginRadius",Math.sqrt(c.p.getOriginX()*c.p.getOriginX()+c.p.getOriginY()+c.p.getOriginY()));
+ instance.addAttValue("particleOriginZ",c.p.getOriginZ());
+ instance.addAttValue("particleStatus",getSimulatorString(c.p.getSimulatorStatus()));
+ instance.addAttValue("split",0);
+
+ return instance;
+ }
+
+}
+
+//convenience class to avoid really long parameter list...
+class InstanceContext
+{
+ public MCParticle p;
+ public ArrayList<SimTrackerHit> hitlist;
+ public HepRepFactory factory;
+ public HepRepInstanceTree instanceTree;
+
+ public InstanceContext(MCParticle part, ArrayList<SimTrackerHit> hits, HepRepFactory fact, HepRepInstanceTree tree) {
+ p = part;
+ hitlist = hits;
+ factory = fact;
+ instanceTree = tree;
+ }
+}
+
+
+//helper class for isCurler()
+class Layer {
+ private int layerNumber;
+ private String detectorName;
+
+ public Layer(int num, String name) {
+ layerNumber=num;
+ detectorName = name;
+ }
+
+ public boolean equals(Object o) {
+ Layer l = (Layer) o;
+ return (l.number()==layerNumber && l.name().equals(detectorName));
+ }
+
+ public int number() {
+ return layerNumber;
+ }
+
+ public String name() {
+ return detectorName;
+ }
+
+}
+
+
+
+
+class TimeCompare<SimTrackerHit> implements Comparator<SimTrackerHit>
+{
+ public int compare(Object a, Object b)
+ {
+ org.lcsim.event.SimTrackerHit hit1 = (org.lcsim.event.SimTrackerHit)a;
+ org.lcsim.event.SimTrackerHit hit2 = (org.lcsim.event.SimTrackerHit)b;
+
+ if (hit1.getTime() < hit2.getTime()) return -1;
+ else if(hit1.getTime() > hit2.getTime()) return 1;
+ else return 0;
+ }
+}
\ No newline at end of file