lcsim/src/org/lcsim/recon/cluster/structural
diff -N ChargedNeutralFragmentSeparator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ChargedNeutralFragmentSeparator.java 3 Oct 2006 00:22:48 -0000 1.1
@@ -0,0 +1,347 @@
+package org.lcsim.recon.cluster.structural;
+
+import java.util.*;
+import hep.physics.vec.*;
+import java.io.IOException;
+
+import org.lcsim.util.decision.DecisionMakerSingle;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.Track;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.Detector;
+
+import hep.aida.ITree;
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogramFactory;
+import hep.aida.ICloud1D;
+
+/**
+ * A class that takes a small cluster or single hit and
+ * tries to decide whether it is charged or neutral.
+ *
+ * @version $Id: ChargedNeutralFragmentSeparator.java,v 1.1 2006/10/03 00:22:48 mcharles Exp $
+ */
+
+public class ChargedNeutralFragmentSeparator extends Driver implements DecisionMakerSingle<Cluster>
+{
+ /** Trivial constructor */
+ public ChargedNeutralFragmentSeparator() {}
+
+ /**
+ * Convenience constructor for the case where we pull down
+ * a new copy of a named list of tracks every event automatically.
+ *
+ * @param trackListName Named list of tracks in the event
+ */
+ public ChargedNeutralFragmentSeparator(String trackListName) {
+ setChargedTracksList(trackListName);
+ }
+
+ /**
+ * Convenience constructor for the case where the user sets the
+ * list of tracks manually.
+ *
+ * @param tracks List of tracks
+ */
+ public ChargedNeutralFragmentSeparator(List<Track> tracks) {
+ setChargedTracksList(tracks);
+ }
+
+ /**
+ * Return true of the object appears to be charged, false if neutral.
+ */
+ public boolean valid (Cluster clus)
+ {
+ // Find the tracks in the event
+ List<Track> tracksToUse = null;
+ if (m_tracksListName != null) {
+ tracksToUse = m_event.get(Track.class, m_tracksListName);
+ } else {
+ tracksToUse = m_tracks;
+ }
+
+ // How many tracks are there in the event?
+ int nTracks = tracksToUse.size();
+ // Solid angle per track
+ double solidAnglePerTrack = Math.PI * 4.0 / ( (double)(nTracks) );
+ // Arbitrary choice of factor... say 5
+ double scaleFactor = 1.5;
+ double coneSolidAngle = scaleFactor * solidAnglePerTrack;
+ double coneOpeningAngle = Math.sqrt(4.0 * coneSolidAngle / Math.PI);
+
+ // Where is the center of the cone?
+ BasicHep3Vector orig = new BasicHep3Vector(0,0,0);
+ BasicHep3Vector clusterPosition = new BasicHep3Vector(clus.getPosition());
+ Hep3Vector clusterDirection = VecOp.sub(clusterPosition, orig);
+ double dotProductThreshold = Math.cos(0.5*coneOpeningAngle);
+
+ // Loop over tracks, checking whether they fall within this cone...
+ List<Track> acceptedTracks = new Vector<Track> ();
+ for (Track tr : tracksToUse) {
+ // It's not sufficient to use the momentum -- need to take the
+ // curvature of the track into account somehow.
+ Hep3Vector swumTrackIntercept = swimTrack(tr);
+ if (swumTrackIntercept != null) {
+ double dotProduct = VecOp.dot(VecOp.unit(clusterDirection), VecOp.unit(swumTrackIntercept));
+ if (dotProduct > dotProductThreshold) {
+ acceptedTracks.add(tr);
+ }
+ }
+ }
+
+ // Use truth information to check whether this cluster is really charged or neutral
+ List<MCParticle> mcList = m_event.get(MCParticle.class, "GenFinalStateParticles");
+ Set<MCParticle> contributingParticles = findMCParticles(clus, mcList);
+ int countChargedContributors = 0;
+ int countNeutralContributors = 0;
+ for (MCParticle part : contributingParticles) {
+ if (Math.abs(part.getCharge()) < 0.5) {
+ countNeutralContributors++;
+ } else {
+ countChargedContributors++;
+ }
+ }
+ String chargeType = new String("error");
+ if (countNeutralContributors>0 && countChargedContributors==0) {
+ chargeType = "neutral";
+ } else if (countNeutralContributors==0 && countChargedContributors>0) {
+ chargeType = "charged";
+ } else if (countNeutralContributors>0 && countChargedContributors>0) {
+ chargeType = "mixed";
+ }
+
+ String printme = new String();
+ printme = "DEBUG: Event has "+tracksToUse.size()+" tracks and scale factor is "
+ +scaleFactor+" so coneSolidAngle="+coneSolidAngle
+ +" and coneOpeningAngle="+coneOpeningAngle
+ +". For a ["+chargeType+"] cluster with "
+ +clus.getCalorimeterHits().size()
+ +" hits, I found "
+ +acceptedTracks.size()+" tracks nearby";
+ if (acceptedTracks.size() >= scaleFactor) {
+ printme += " => charged";
+ } else {
+ printme += " => neutral";
+ }
+ System.out.println(printme);
+
+ boolean recoIsCharged = acceptedTracks.size() >= scaleFactor;
+ boolean trueIsCharged = (countChargedContributors>0 && countNeutralContributors==0);
+ if (recoIsCharged == trueIsCharged) {
+ m_fragPerformance.fill(1);
+ } else {
+ m_fragPerformance.fill(0);
+ }
+ if (trueIsCharged) {
+ m_fragPerformanceCharged.fill( (double)(acceptedTracks.size()) / scaleFactor );
+ } else {
+ m_fragPerformanceNeutral.fill( (double)(acceptedTracks.size()) / scaleFactor );
+ }
+
+ return (acceptedTracks.size() >= scaleFactor);
+ }
+
+ public void setChargedTracksList(String listName) {
+ m_tracksListName = listName;
+ m_tracks = null;
+ }
+ public void setChargedTracksList(List<Track> tracks) {
+ m_tracksListName = null;
+ m_tracks = tracks;
+ }
+ public void process(EventHeader event) {
+ m_event = event;
+ initGeometry(event);
+ Detector det = event.getDetector();
+ double[] zero = {0, 0, 0};
+ double[] fieldStrength = det.getFieldMap().getField(zero);
+ m_swimmer = new HelixSwimmer(fieldStrength[2]);
+ }
+
+ protected boolean m_init = false;
+ protected HelixSwimmer m_swimmer = null;
+ protected EventHeader m_event = null;
+ protected String m_tracksListName = null;
+ protected List<Track> m_tracks = null;
+
+ protected double m_ECAL_barrel_zmin;
+ protected double m_ECAL_barrel_zmax;
+ protected double m_ECAL_barrel_r;
+ protected double m_ECAL_endcap_z;
+ protected double m_ECAL_endcap_rmin;
+ protected double m_ECAL_endcap_rmax;
+
+ // FOR DEBUG
+
+ ITree m_tree;
+ IHistogramFactory m_histoFactory = null;
+ ICloud1D m_fragPerformance;
+ ICloud1D m_fragPerformanceNeutral;
+ ICloud1D m_fragPerformanceCharged;
+ public void initPlots()
+ {
+ IAnalysisFactory af = IAnalysisFactory.create();
+ try {
+ m_tree = af.createTreeFactory().create("charged-neutral.aida","xml",false,true);
+ m_histoFactory = af.createHistogramFactory(m_tree);
+ m_fragPerformance = m_histoFactory.createCloud1D("fragPerformance");
+ m_fragPerformanceNeutral = m_histoFactory.createCloud1D("fragPerformanceNeutral");
+ m_fragPerformanceCharged = m_histoFactory.createCloud1D("fragPerformanceCharged");
+ } catch (IOException ioe1) {
+ ioe1.printStackTrace();
+ }
+ }
+
+ public void suspend() {
+ try {
+ m_tree.commit();
+ } catch(IOException ioe1) {
+ ioe1.printStackTrace();
+ }
+ super.suspend();
+ }
+
+ /**
+ * Internal utility routine, belongs in another class
+ */
+ protected Set<MCParticle> findMCParticles(Cluster clus, List<MCParticle> mcList)
+ {
+ Set clusterParticles = new HashSet<MCParticle>();
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ Set<MCParticle> hitParticles = findMCParticles(hit, mcList);
+ clusterParticles.addAll(hitParticles);
+ }
+ return clusterParticles;
+ }
+
+ /**
+ * Internal utility routine, belongs in another class
+ */
+ protected Set<MCParticle> findMCParticles(CalorimeterHit hit, List<MCParticle> mcList)
+ {
+ if ( ! (hit instanceof SimCalorimeterHit) ) {
+ throw new AssertionError("Non-simulated hit!");
+ } else {
+ SimCalorimeterHit simHit = (SimCalorimeterHit) (hit);
+ Set<MCParticle> contributingParticlesFromList = new HashSet<MCParticle>();
+ int nContributingParticles = simHit.getMCParticleCount();
+ for (int i=0; i<nContributingParticles; i++) {
+ MCParticle part = simHit.getMCParticle(i);
+ List<MCParticle> parentsInList = findParentsInList(part, mcList);
+ contributingParticlesFromList.addAll(parentsInList);
+ }
+ return contributingParticlesFromList;
+ }
+ }
+
+ /**
+ * Internal utility routine, belongs in another class
+ */
+ protected List<MCParticle> findParentsInList(MCParticle part, List<MCParticle> mcList)
+ {
+ List<MCParticle> outputList = new Vector<MCParticle>();
+ if (mcList.contains(part)) {
+ // Already in there
+ outputList.add(part);
+ } else {
+ // Not in there -- recurse up through parents
+ List<MCParticle> parents = part.getParents();
+ if (parents.size()==0) {
+ // Ran out of options -- add nothing and return below
+ } else {
+ for (MCParticle parent : parents) {
+ List<MCParticle> ancestorsInList = findParentsInList(parent, mcList);
+ outputList.addAll(ancestorsInList);
+ }
+ }
+ }
+ return outputList;
+ }
+
+ protected void initGeometry(EventHeader event)
+ {
+ if (!m_init) {
+ Detector det = event.getDetector();
+ CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
+ CylindricalCalorimeter eme = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap"));
+ m_ECAL_barrel_zmin = emb.getZMin();
+ m_ECAL_barrel_zmax = emb.getZMax();
+ m_ECAL_barrel_r = emb.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_z = eme.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_rmin = eme.getInnerRadius();
+ m_ECAL_endcap_rmax = eme.getOuterRadius();
+ m_init = true;
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_zmin="+m_ECAL_barrel_zmin);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_zmax="+m_ECAL_barrel_zmax);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_r="+m_ECAL_barrel_r);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_z="+m_ECAL_endcap_z);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_rmin="+m_ECAL_endcap_rmin);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_rmax="+m_ECAL_endcap_rmax);
+ }
+ }
+
+ // Swim a track
+ protected Hep3Vector swimTrack(Track tr) {
+ double alpha = Double.NaN;
+ m_swimmer.setTrack(tr);
+ // Try swimming to the barrel:
+ double alphaBarrel = swimToBarrel(m_swimmer);
+ boolean validBarrel = false;
+ // Try swimming to the endcap:
+ double alphaEndcap = swimToEndcap(m_swimmer);
+ boolean validEndcap = false;
+ // Fixme: Here we should check that the track really does go all the
+ // way to the ECAL instead of stopping/decaying/interacting earlier.
+ // This used to be done with the checkDecayPoint() method in
+ // contrib.uiowa.structural.SwimToECAL
+ if (isValidBarrelIntercept(m_swimmer, alphaBarrel)) {
+ alpha = alphaBarrel;
+ validBarrel = true;
+ } else if (isValidEndcapIntercept(m_swimmer, alphaEndcap)) {
+ alpha = alphaEndcap;
+ validEndcap = true;
+ }
+
+ if ( ! Double.isNaN(alpha) ) {
+ if (validEndcap || validBarrel) {
+ if ( !(validEndcap && validBarrel) ) {
+ // OK
+ Hep3Vector trackPoint = m_swimmer.getPointAtDistance(alpha);
+ return trackPoint;
+ }
+ }
+ }
+
+ // Unable to extrapolate
+ return null;
+ }
+
+ protected double swimToBarrel(HelixSwimmer swimmer) {
+ // Look for a hit in the first layer of the ECAL barrel
+ return swimmer.getDistanceToRadius(m_ECAL_barrel_r);
+ }
+ protected double swimToEndcap(HelixSwimmer swimmer) {
+ // Look for a hit in the first layer of the ECAL endcap
+ return swimmer.getDistanceToZ(m_ECAL_endcap_z);
+ }
+ protected boolean isValidBarrelIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have -m_ECAL_barrel_z <= z <= +m_ECAL_barrel_z
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double z = intercept.z();
+ boolean zInRange = (z >= m_ECAL_barrel_zmin && z <= m_ECAL_barrel_zmax);
+ return zInRange;
+ }
+ protected boolean isValidEndcapIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have m_ECAL_endcap_rmin <= r <= m_ECAL_endcap_rmax
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double r = Math.sqrt(intercept.x()*intercept.x() + intercept.y()*intercept.y());
+ boolean rInRange = (r >= m_ECAL_endcap_rmin && r <= m_ECAL_endcap_rmax);
+ return rInRange;
+ }
+}
