diff -u -r1.7 -r1.8
--- GBLOutput.java	25 Aug 2013 21:42:19 -0000	1.7
+++ GBLOutput.java	26 Aug 2013 21:40:03 -0000	1.8
@@ -46,7 +46,8 @@

     private MaterialSupervisor _materialmanager;
     private MultipleScattering _scattering;
     private Detector detector;

-    private HPSTransformations _hpstrans = new HPSTransformations();                

+    private HPSTransformations _hpstrans = new HPSTransformations();
+    private double _beamEnergy = 2.2; //GeV

@@ -123,9 +124,34 @@

         // Use the sim tracker hits to find the truth particle
         // this would give me the particle kinematics responsible for the track 
         

+        List<MCParticle> mcp_pair = getAprimeDecayProducts(mcParticles);
+        if(_debug>0) {
+            double invMassTruth = Math.sqrt( Math.pow(mcp_pair.get(0).getEnergy()+mcp_pair.get(1).getEnergy(),2) - VecOp.add(mcp_pair.get(0).getMomentum(), mcp_pair.get(1).getMomentum()).magnitudeSquared() );
+            double invMassTruthTrks = getInvMassTracks(getHTF(mcp_pair.get(0)),getHTF(mcp_pair.get(1)));
+            System.out.printf("%s: invM = %f\n",this.getClass().getSimpleName(),invMassTruth);
+            System.out.printf("%s: invMTracks = %f\n",this.getClass().getSimpleName(),invMassTruthTrks);
+        }

         
         MCParticle mcp = getMatchedTruthParticle(trk);

-

+        
+        // cross-check
+        if(!mcp_pair.contains(mcp)) {
+            boolean hasBeamElectronParent = false;
+            for(MCParticle parent : mcp.getParents()) {
+                if(parent.getGeneratorStatus()!=MCParticle.FINAL_STATE && parent.getPDGID()==11 && parent.getMomentum().y()==0.0 && Math.abs(parent.getMomentum().magnitude() - _beamEnergy)<0.01) {
+                    hasBeamElectronParent = true;
+                }
+            }
+            if(!hasBeamElectronParent) {
+                System.out.printf("%s: the matched MC particle is not an A' daughter and not a the recoil electrons!?\n",this.getClass().getSimpleName());
+                System.out.printf("%s: %s %d p %s org %s\n",this.getClass().getSimpleName(),mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString());
+                printMCParticles(mcParticles);
+                System.exit(1);
+            } else {
+                if(_debug>0) System.out.printf("%s: the matched MC particle is the recoil electron\n",this.getClass().getSimpleName());
+            }
+        }
+        

         List<SimTrackerHit> simHits = new ArrayList<SimTrackerHit>();
         for (SimTrackerHit hit : simTrackerHits) {
             if(hit.getMCParticle()==mcp) simHits.add(hit);

@@ -214,7 +240,7 @@

                 }
                 
                 Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z()));

-

+                

                 if(_debug>0) {
                     System.out.printf("trkposTruth at intercept %s\n",trkposTruth.toString());
                 }

@@ -247,7 +273,13 @@

                         System.out.printf("\n");
                     }
                     //System.exit(1);

-                } 

+                }
+                if(_debug>0) {
+                    double s_truthSimHit = HelixUtils.PathToXPlane(htfTruth, simHit.getPositionVec().z(), 0, 0).get(0);
+                    Hep3Vector trkposTruthSimHit = HelixUtils.PointOnHelix(htfTruth, s_truthSimHit);
+                    Hep3Vector resTruthSimHit = VecOp.sub(this._hpstrans.transformVectorToTracking(simHit.getPositionVec()),trkposTruthSimHit);
+                    System.out.printf("TruthSimHit residual %s for layer %d\n",resTruthSimHit.toString(),strip.layer());
+                }

                 
                 //path length to intercept
                 double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);   

@@ -649,9 +681,74 @@

         return chi2.e(0, 0);
     }
 

+    
+    
+    private List<MCParticle> getAprimeDecayProducts(List<MCParticle> mcParticles) {
+        List<MCParticle> pair = new ArrayList<MCParticle>();
+        for(MCParticle mcp : mcParticles) {
+            if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
+            boolean hasAprimeParent = false;
+            for(MCParticle parent : mcp.getParents()) {
+                if(Math.abs(parent.getPDGID())==622) hasAprimeParent = true;
+            }
+            if(hasAprimeParent)  pair.add(mcp);
+        }
+        if(pair.size()!=2) {
+            System.out.printf("%s: ERROR this event has %d mcp with 622 as parent!!??  \n",this.getClass().getSimpleName(),pair.size());
+            this.printMCParticles(mcParticles);
+            System.exit(1);
+        }
+        if( Math.abs(pair.get(0).getPDGID()) != 11 || Math.abs(pair.get(1).getPDGID()) != 11 ) {
+            System.out.printf("%s: ERROR decay products are not e+e-? \n",this.getClass().getSimpleName());
+            this.printMCParticles(mcParticles);
+            System.exit(1);
+        }
+        if(pair.get(0).getPDGID()*pair.get(1).getPDGID() > 0) {
+            System.out.printf("%s: ERROR decay products have the same sign? \n",this.getClass().getSimpleName());
+            this.printMCParticles(mcParticles);
+            System.exit(1);
+        }
+        return pair;
+        
+    }
+
+    private void printMCParticles(List<MCParticle> mcParticles) {      
+        System.out.printf("%s: printMCParticles \n",this.getClass().getSimpleName());
+        System.out.printf("%s: %d mc particles \n",this.getClass().getSimpleName(),mcParticles.size());
+        for(MCParticle mcp : mcParticles) {
+            if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
+            System.out.printf("\n%s: (%s) %d  p %s org %s  %s \n",this.getClass().getSimpleName(),
+                                mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString(),
+                                mcp.getParents().size()>0?"parents:":"");
+            for(MCParticle parent : mcp.getParents()) {
+                System.out.printf("%s:       (%s) %d  p %s org %s %s \n",this.getClass().getSimpleName(),
+                                parent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",parent.getPDGID(),parent.getMomentum().toString(),parent.getOrigin().toString(),
+                                parent.getParents().size()>0?"parents:":"");
+                for(MCParticle grparent : parent.getParents()) {
+                    System.out.printf("%s:            (%s) %d  p %s org %s  %s \n",this.getClass().getSimpleName(),
+                                    grparent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",grparent.getPDGID(),grparent.getMomentum().toString(),grparent.getOrigin().toString(),
+                                    grparent.getParents().size()>0?"parents:":"");
+                }
+            
+            }
+        }
+        return;
+     }

-   

+    private double getInvMassTracks(HelicalTrackFit htf1, HelicalTrackFit htf2) {
+        double p1 = htf1.p(this._B.magnitude());
+        double p2 = htf2.p(this._B.magnitude());
+        Hep3Vector p1vec = VecOp.mult(p1, HelixUtils.Direction(htf1, 0));
+        Hep3Vector p2vec = VecOp.mult(p2, HelixUtils.Direction(htf2, 0));
+        double me = 0.000510998910;
+        double E1 = Math.sqrt(p1*p1 + me*me);
+        double E2 = Math.sqrt(p2*p2 + me*me);
+        //System.out.printf("p1 %f %s E1 %f\n",p1,p1vec.toString(),E1);
+        //System.out.printf("p2 %f %s E2 %f\n",p2,p2vec.toString(),E2);
+        return Math.sqrt( Math.pow(E1+E2,2) - VecOp.add(p1vec, p2vec).magnitudeSquared() );
+    }

+    

     public class PerigeeParams {
         private BasicMatrix _params;