Author: [log in to unmask]
Date: Fri May  1 17:11:04 2015
New Revision: 2876

Log:
add track truth mc function. add debug option.

Modified:
    java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java

Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
 =============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java	(original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java	Fri May  1 17:11:04 2015
@@ -192,7 +192,7 @@
      * @return position at ECAL
      */
     public static Hep3Vector getTrackPositionAtEcal(Track track) {
-        return extrapolateTrack(track, BeamlineConstants.ECAL_FACE_TESTRUN);
+        return extrapolateTrack(track, BeamlineConstants.ECAL_FACE);
     }
 
     /**
@@ -663,8 +663,18 @@
      * @return helix object based on the MC particle
      */
     public static HelicalTrackFit getHTF(MCParticle mcp, double Bz) {
+        boolean debug = true;
+        if(debug) {
+            System.out.printf("getHTF\n");
+            System.out.printf("mcp org %s mc p %s\n",mcp.getOrigin().toString(),mcp.getMomentum().toString());
+        }
         Hep3Vector org = CoordinateTransformations.transformVectorToTracking(mcp.getOrigin());
         Hep3Vector p = CoordinateTransformations.transformVectorToTracking(mcp.getMomentum());
+
+        if(debug) {
+            System.out.printf("mcp org %s mc p %s (trans)\n",org.toString(),p.toString());
+        }
+
         // Move to x=0 if needed
         double targetX = BeamlineConstants.DIPOLE_EDGELOW_TESTRUN;
         if (org.x() < targetX) {
@@ -681,7 +691,9 @@
             // old.toString(),p.toString(),org.toString());
         }
 
-        // System.out.printf("outside org %s p %s \n",p.toString(),org.toString());
+        if(debug) {
+            System.out.printf("mcp org %s mc p %s (trans2)\n",org.toString(),p.toString());
+        }
 
         HelixParamCalculator helixParamCalculator = new HelixParamCalculator(p, org, -1 * ((int) mcp.getCharge()), Bz);
         double par[] = new double[5];
@@ -691,8 +703,8 @@
         par[HelicalTrackFit.curvatureIndex] = 1.0 / helixParamCalculator.getRadius();
         par[HelicalTrackFit.z0Index] = helixParamCalculator.getZ0();
         HelicalTrackFit htf = getHTF(par);
-        // System.out.printf("d0 %f z0 %f R %f phi %f lambda %s\n",
-        // htf.dca(),htf.z0(),htf.R(),htf.phi0(),htf.slope() );
+         System.out.printf("d0 %f z0 %f R %f phi %f lambda %s\n",
+         htf.dca(),htf.z0(),htf.R(),htf.phi0(),htf.slope() );
         return htf;
     }
 
@@ -735,5 +747,59 @@
         // System.out.printf("%s: straight line track: x0=%f,y0=%f,z0=%f dz/dx=%f dydx=%f targetY=%f targetZ=%f \n",this.getClass().getSimpleName(),slt1.x0(),slt1.y0(),slt1.z0(),slt1.dzdx(),slt1.dydx(),slt1.TargetYZ()[0],slt1.TargetYZ()[1]);
         return slt1;
     }
+    
+    
+    
+    public static MCParticle getMatchedTruthParticle(Track track) {
+        boolean debug = false;
+        
+        Map<MCParticle,Integer> particlesOnTrack = new HashMap<MCParticle,Integer>();
+        
+        if(debug) System.out.printf("getMatchedTruthParticle: getmatched mc particle from %d tracker hits on the track \n",track.getTrackerHits().size());
+        
+        
+        for(TrackerHit hit : track.getTrackerHits()) {
+            List<MCParticle> mcps = ((HelicalTrackHit)hit).getMCParticles();
+            if(mcps==null) {
+                System.out.printf("getMatchedTruthParticle: warning, this hit (layer %d pos=%s) has no mc particles.\n",((HelicalTrackHit)hit).Layer(),((HelicalTrackHit)hit).getCorrectedPosition().toString());
+            } 
+            else {
+                if( debug ) System.out.printf("getMatchedTruthParticle: this hit (layer %d pos=%s) has %d mc particles.\n",((HelicalTrackHit)hit).Layer(),((HelicalTrackHit)hit).getCorrectedPosition().toString(),mcps.size());
+                for(MCParticle mcp : mcps) {
+                    if( !particlesOnTrack.containsKey(mcp) ) {
+                        particlesOnTrack.put(mcp, 0);
+                    }
+                    int c = particlesOnTrack.get(mcp);
+                    particlesOnTrack.put(mcp, c+1);
+                }
+            }
+        }
+        if(debug) {
+            System.out.printf("Track p=[ %f, %f, %f] \n",track.getTrackStates().get(0).getMomentum()[0],track.getTrackStates().get(0).getMomentum()[1],track.getTrackStates().get(0).getMomentum()[1]);
+            System.out.printf("Found %d particles\n",particlesOnTrack.size());
+            for(Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
+                System.out.printf("%d hits assigned to %d p=%s \n",entry.getValue(),entry.getKey().getPDGID(),entry.getKey().getMomentum().toString());
+            }
+        }
+        Map.Entry<MCParticle,Integer> maxEntry = null;
+        for(Map.Entry<MCParticle,Integer> entry : particlesOnTrack.entrySet()) {
+            if ( maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0 ) maxEntry = entry;
+            //if ( maxEntry != null ) {
+            //    if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue;
+            //}
+            //maxEntry = entry;
+        }
+        if(debug) {
+            if (maxEntry != null ) {
+                System.out.printf("Matched particle with pdgId=%d and mom %s to track with charge %d and momentum [%f %f %f]\n",
+                        maxEntry.getKey().getPDGID(),maxEntry.getKey().getMomentum().toString(),
+                        track.getCharge(),track.getTrackStates().get(0).getMomentum()[0],track.getTrackStates().get(0).getMomentum()[1],track.getTrackStates().get(0).getMomentum()[2]);
+            } else {
+                System.out.printf("No truth particle found on this track\n");
+            }
+        }
+        return maxEntry == null ? null : maxEntry.getKey();
+    }
+    
 
 }