lcsim/src/org/lcsim/contrib/timb/mc/fast/tracking
diff -N MCFastTracking.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCFastTracking.java 26 May 2006 07:21:59 -0000 1.1
@@ -0,0 +1,178 @@
+package org.lcsim.mc.fast.tracking;
+
+import hep.physics.particle.Particle;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.Iterator;
+import java.util.List;
+
+import org.lcsim.conditions.ConditionsEvent;
+import org.lcsim.conditions.ConditionsListener;
+import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.util.Driver;
+import org.lcsim.particle.ParticleUtilities;
+
+
+/**
+ * Fast Monte Carlo tracking simulator
+ */
+public class MCFastTracking extends Driver implements ConditionsListener
+{
+ private List<Particle> bquark_list = new ArrayList<Particle>();
+ private List<Particle> anti_bquark_list = new ArrayList<Particle>();
+ private Particle cquark;
+ private Particle anti_cquark;
+ private TrackResolutionTables parm;
+ private SimpleTables SmTbl;
+ private boolean beamSpotConstraint;
+ private boolean simple;
+ private final static double[] IP = { 0, 0, 0 };
+ private static int nprint = 0 ;
+ private static final int nprint_max = -2 ;
+
+ public MCFastTracking()
+ {
+ this(false);
+ }
+
+ public MCFastTracking(boolean beamSpotConstraint)
+ {
+ this.beamSpotConstraint = beamSpotConstraint;
+ }
+
+ public MCFastTracking(boolean beamSpotConstraint, boolean simple)
+ {
+ this.beamSpotConstraint = beamSpotConstraint;
+ this.simple = simple;
+ }
+
+ public void setBeamSpotConstraint(boolean beamSpotConstraint)
+ {
+ this.beamSpotConstraint = beamSpotConstraint;
+ if (parm != null)
+ {
+ ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
+ parm = setTrackResolutionTables(conditions,beamSpotConstraint);
+ }
+ }
+
+ public boolean isBeamSpotConstraint()
+ {
+ return this.beamSpotConstraint;
+ }
+ private TrackResolutionTables setTrackResolutionTables(ConditionsSet conditions,boolean beamSpotConstraint)
+ {
+ try
+ {
+ return new TrackResolutionTables(conditions, beamSpotConstraint);
+ }
+ catch (IOException x)
+ {
+ throw new RuntimeException("Error reading track resolution tables",x);
+ }
+ }
+
+ protected void process(EventHeader event)
+ {
+ if (parm == null)
+ {
+ ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
+ conditions.addConditionsListener(this);
+ parm = setTrackResolutionTables(conditions, beamSpotConstraint);
+ }
+
+ if (SmTbl == null)
+ {
+ ConditionsSet simpleconditions = getConditionsManager().getConditions("SimpleTrack");
+ simpleconditions.addConditionsListener(this);
+ SmTbl = new SimpleTables(simpleconditions);
+ }
+
+ double bField = event.getDetector().getFieldMap().getField(IP)[2];
+ boolean hist = getHistogramLevel() > 0;
+
+ List<Track> trackList = new ArrayList<Track>();
+ bquark_list.clear();
+ anti_bquark_list.clear();
+ cquark=null;
+ anti_cquark=null;
+ for (Iterator i = event.getMCParticles().iterator(); i.hasNext();)
+ {
+ Particle p = (Particle) i.next();
+ if(Math.abs(p.getPDGID()) == 5) {
+ boolean addp;
+
+ if(p.getParents().size() == 0) addp = false;
+ else addp=Math.abs(((Particle)p.getParents().get(0)).getPDGID()) != 5;
+
+ if(p.getPDGID() == 5 && addp) bquark_list.add(p);
+ if(p.getPDGID() == -5 && addp) anti_bquark_list.add(p);
+ }
+ if(p.getPDGID() == 4) cquark=p;
+ if(p.getPDGID() == -4) anti_cquark=p;
+
+ // filter for FINAL_STATE
+ if (p.getGeneratorStatus() != Particle.FINAL_STATE)
+ {
+ continue;
+ }
+ double pCharge = p.getCharge();
+ if (pCharge == 0 || pCharge == Double.NaN || pCharge == Double.NEGATIVE_INFINITY || pCharge == Double.POSITIVE_INFINITY)
+ {
+ continue;
+ }
+
+ double[] momentum = p.getMomentum().v();
+ double pt2 = (momentum[0] * momentum[0]) + (momentum[1] * momentum[1]);
+ double pt = Math.sqrt(pt2);
+ double ptot = Math.sqrt(pt2 + (momentum[2] * momentum[2]));
+ double cosTheta = momentum[2] / ptot;
+
+ // within acceptance
+ if (pt < parm.getPtMin())
+ {
+ continue;
+ }
+ if (Math.abs(cosTheta) > parm.getPolarOuter())
+ {
+ continue;
+ }
+
+ trackList.add(new ReconTrack(bField, parm, SmTbl, getRandom(), p, hist, simple));
+ }
+ nprint++;
+ if(nprint <= nprint_max) {
+
+ for(Particle bquark : bquark_list) {
+ System.out.println(" MCFastTracking nprint= "+nprint+" bquark.getPDGID= "+bquark.getPDGID()+" bquark.getGeneratorStatus= "+bquark.getGeneratorStatus()
+ +" ((Particle)bquark.getParents().get(0)).getPDGID= "+((Particle)bquark.getParents().get(0)).getPDGID());
+ ParticleUtilities.dumpParticleHierarchy(bquark);
+ }
+ for(Particle anti_bquark : anti_bquark_list) {
+ System.out.println(" MCFastTracking nprint= "+nprint+" anti_bquark.getPDGID= "+anti_bquark.getPDGID()+" anti_bquark.getGeneratorStatus= "+anti_bquark.getGeneratorStatus()
+ +" ((Particle)anti_bquark.getParents().get(0)).getPDGID= "+((Particle)anti_bquark.getParents().get(0)).getPDGID());
+ ParticleUtilities.dumpParticleHierarchy(anti_bquark);
+ }
+ if(cquark != null && bquark_list.size() == 0) {
+ System.out.println(" MCFastTracking nprint= "+nprint+" cquark.getPDGID= "+cquark.getPDGID()+" cquark.getGeneratorStatus= "+cquark.getGeneratorStatus());
+ ParticleUtilities.dumpParticleHierarchy(cquark);
+ }
+ if(anti_cquark != null && anti_bquark_list.size() == 0) {
+ System.out.println(" MCFastTracking nprint= "+nprint+" anti_cquark.getPDGID= "+anti_cquark.getPDGID()+" anti_cquark.getGeneratorStatus= "+anti_cquark.getGeneratorStatus());
+ ParticleUtilities.dumpParticleHierarchy(anti_cquark);
+ }
+ }
+
+ event.put(EventHeader.TRACKS, trackList, Track.class, 0);
+ }
+
+ public void conditionsChanged(ConditionsEvent event)
+ {
+ ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
+ ConditionsSet simpleconditions = getConditionsManager().getConditions("SimpleTrack");
+ parm = setTrackResolutionTables(conditions, beamSpotConstraint);
+ SmTbl = new SimpleTables(simpleconditions);
+ }
+}