hps-java/src/main/java/org/lcsim/hps/users/phansson

diff -u -r1.8 -r1.9
--- RunMPAlignment.java	12 Oct 2012 05:59:36 -0000	1.8
+++ RunMPAlignment.java	23 Oct 2012 01:31:17 -0000	1.9
@@ -25,8 +25,6 @@

 
     private AIDA aida = AIDA.defaultInstance();
     String[] detNames = {"Tracker"};

-    Integer _minLayers = 8;
-    Integer[] nlayers = {8};

     int nevt = 0;
     double[] beamsize = {0.001, 0.02, 0.02};
     String _config = "";

@@ -43,16 +41,17 @@

 //  this is set by the _config variable (detType in HeavyPhotonDriver)
     int flipSign = 1;
     private boolean _debug = false;

-    double _resLayer1MinY;
-    

+    private String _type = "LOCAL"; //GLOBAL OR LOCAL RESIDUALS

-   String simTrackerHitCollectionName = "TrackerHits";
-   

+    private String simTrackerHitCollectionName = "TrackerHits";
+    private String outputMilleFile = "alignMP.txt";

     
     public void setDebug(boolean v) {
         this._debug = v;
     }

-    

+    public void setMilleFile(String filename) {
+        outputMilleFile = filename;
+    }

     public void setOutputPlotFileName(String filename) {
         outputPlotFileName = filename;
     }

@@ -61,60 +60,35 @@

         hideFrame = hide;
     }
     

-    public RunMPAlignment() {
-        nlayers[0] = 10;
-        _minLayers = 8;
-////        if (_config.contains("Test"))
-////            flipSign = -1;
- 
-
-    

+    public void setResidualLimitFileName(String fileName) {
+        this._resLimitFileName = fileName;

     }

-    
-    public RunMPAlignment(int trackerLayers, int mintrkLayers, String config) {
-        nlayers[0] = trackerLayers;
-        _minLayers = mintrkLayers;
-        _config = config;
-        if (_config.contains("Test"))
-            flipSign = -1;
-        

+    
+    public RunMPAlignment() {

     }
 
     
     @Override
     public void detectorChanged(Detector detector) {
         

-        //ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug,hideFrame);
-        String outputMilleTxtFile = "alignMP.txt";
-        Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
-        Hep3Vector _bfield = new BasicHep3Vector(0,0,detector.getFieldMap().getField(IP).y());
-
-        ap = new MPAlignmentParameters(outputMilleTxtFile,_bfield,_debug,hideFrame);
-               
-        //Initialize the res limits
-        for(int i=1;i<=10;++i) {
-            for(int j=0;j<3;++j) {
-                double xmin = -50;
-                double xmax = 50;
-                for(int side=0;side<2;++side) {
-                   ap.setResLimits(side,i,j, xmin, xmax);        
-                }
-            }
-        }

+        ap = new MPAlignmentParameters(this.outputMilleFile);
+        ap._DEBUG = _debug;
+        ap.setType(this._type);
+        ap.setHideFrame(hideFrame);
+        double bfield = detector.getFieldMap().getField(new BasicHep3Vector(0., 0., 1.)).y();
+        System.out.printf("%s: B-field in z %.3f\n",this.getClass().getSimpleName(),bfield);
+        ap.setBfield(bfield);

         loadResidualLimits();
         

-        
-        

     }
     
     
     
     @Override

-    public void process(
-            EventHeader event) {
-

+    public void process(EventHeader event) {

+        

         List<Track> tracklist = null;
         if(event.hasCollection(Track.class,"MatchedTracks")) {        
             tracklist = event.get(Track.class, "MatchedTracks");

@@ -123,33 +97,6 @@

              }
         }
         

-        List<SimTrackerHit> simTrackerHits = new ArrayList<SimTrackerHit>();
-        if(event.hasCollection(SimTrackerHit.class, simTrackerHitCollectionName)){
-            simTrackerHits.addAll(event.get(SimTrackerHit.class, simTrackerHitCollectionName));
-            if(_debug){
-                System.out.println(this.getClass().getSimpleName() + ": Event: " + event.getEventNumber() + " Number of SimTrackerHits:  " + simTrackerHits.size());
-                System.out.print(this.getClass().getSimpleName() + ": MC Particles: ");
-                for(MCParticle mcParticle : event.getMCParticles()){
-                    System.out.print(mcParticle.getPDGID() + " ");
-                }
-                System.out.print("\n");
-            }
-        }
-        List<HelicalTrackStrip> strips = null;
-        if(event.hasCollection(HelicalTrackStrip.class, "HelicalTrackStrips")) {
-            strips = event.get(HelicalTrackStrip.class, "HelicalTrackStrips");
-            if(_debug) System.out.println(this.getClass().getSimpleName() + ": Event has " + strips.size() + " HelicalTrackStrips");
-        }
-
-        List<SiTrackerHit> trackerHits = null;
-                
-        if(event.hasCollection(SiTrackerHit.class, "StripClusterer_SiTrackerHitStrip1D")) {
-            trackerHits = event.get(SiTrackerHit.class, "StripClusterer_SiTrackerHitStrip1D");
-            if(_debug) System.out.println(this.getClass().getSimpleName() + ": Event has " + trackerHits.size() + " SiTrackerHit");
-        }
-
-        
-       

         for (Track trk : tracklist) {
             
             //if(trk.getCharge()>0) continue;

@@ -159,23 +106,28 @@

             
             if(hitsOnBothSides(trk)) continue;
             

-            totalTracksProcessed++;

+            double Px =  trk.getTrackStates().get(0).getMomentum()[0];
+            if(Px < 0.1) {
+                System.out.printf("%s: Trk p = [%.3f,%.3f,%.3f] is low skip!?\n",this.getClass().getSimpleName(),
+                        trk.getTrackStates().get(0).getMomentum()[0],trk.getTrackStates().get(0).getMomentum()[1],trk.getTrackStates().get(0).getMomentum()[2]);
+                continue;
+            }

+            totalTracksProcessed++;

-            ap.PrintResidualsAndDerivatives(trk,totalTracks,"LOCAL",simTrackerHits);

+            ap.PrintResidualsAndDerivatives(trk, totalTracks);

             
             
             
 
         }

-
-        ap.CalculateResidualSimHits(strips,trackerHits,simTrackerHits);

         
         
         
         
     }
 

+    @Override

     public void endOfData() {
         ap.updatePlots();
         try {

@@ -185,11 +137,12 @@

         } catch (IOException ex) {
             Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE, null, ex);
         }

-        if (outputPlotFileName != "")
-        try {
-            aida.saveAs(outputPlotFileName);
-        } catch (IOException ex) {
-            Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);

+        if (!"".equals(outputPlotFileName)) {
+            try {
+                aida.saveAs(outputPlotFileName);
+            } catch (IOException ex) {
+                Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
+            }

         }
         
     }

@@ -220,49 +173,58 @@

     
     private void loadResidualLimits() {
         

-        FileReader fReader = null;
-        BufferedReader bReader = null;
-        try {
-          fReader = new FileReader(this._resLimitFileName);
-          bReader = new BufferedReader(fReader);
-          
-          String line;
-          while( (line = bReader.readLine()) != null) {
-              if (line.contains("#")) continue;
-              String[] vec = line.split("\\s+");
-              if(vec.length!=5) {
-                  System.out.println(this.getClass().getSimpleName() + ": Error: residual limits line has wrong format -> " + line);
-                  System.exit(1);
-              }
-              try {
-                int side = Integer.parseInt(vec[0]);
-                int layer = Integer.parseInt(vec[1]);
-                int direction = Integer.parseInt(vec[2]);
-                double min = Double.parseDouble(vec[3]);
-                double max = Double.parseDouble(vec[4]);
-                ap.setResLimits(side, layer, direction, min, max);
-              } catch(NumberFormatException e) {
-                  Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE,null,e);
-              }
-          }
-          bReader.close();
-          fReader.close();
-        }
-        catch (FileNotFoundException ex) {
-            Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE, null, ex);
-        } catch (IOException e) {
-            Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE,null,e);
-        } 

+        //Initialize the res limits
+        for(int i=1;i<=10;++i) {
+            for(int j=0;j<3;++j) {
+                double xmin = -50;
+                double xmax = 50;
+                for(int side=0;side<2;++side) {
+                   ap.setResLimits(side,i,j, xmin, xmax);        
+                }
+            }
+        }

+        if(!"".equals(this._resLimitFileName)) {
+            FileReader fReader;
+            BufferedReader bReader;
+            try {
+                fReader = new FileReader(this._resLimitFileName);
+                bReader = new BufferedReader(fReader);
+
+                String line;
+                while( (line = bReader.readLine()) != null) {
+                    if (line.contains("#")) continue;
+                    String[] vec = line.split("\\s+");
+                    if(vec.length!=5) {
+                        System.out.println(this.getClass().getSimpleName() + ": Error: residual limits line has wrong format -> " + line);
+                        System.exit(1);
+                    }
+                    try {
+                        int side = Integer.parseInt(vec[0]);
+                        int layer = Integer.parseInt(vec[1]);
+                        int direction = Integer.parseInt(vec[2]);
+                        double min = Double.parseDouble(vec[3]);
+                        double max = Double.parseDouble(vec[4]);
+                        ap.setResLimits(side, layer, direction, min, max);
+                    } catch(NumberFormatException e) {
+                        Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE,null,e);
+                    }
+                }
+                bReader.close();
+                fReader.close();
+            }
+            catch (FileNotFoundException ex) {
+                Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE, null, ex);
+            } catch (IOException e) {
+                Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE,null,e);
+            } 
+        }

-        ap.setResLimits(0,1,0, -0.3, 0.3);

+        

 
     }
     
     

-    public void setResidualLimitFileName(String fileName) {
-        this._resLimitFileName = fileName;
-    }

hps-java/src/main/java/org/lcsim/hps/users/phansson

diff -u -r1.13 -r1.14
--- MPAlignmentParameters.java	17 Oct 2012 23:00:38 -0000	1.13
+++ MPAlignmentParameters.java	23 Oct 2012 01:31:17 -0000	1.14
@@ -2,34 +2,29 @@

 
 import hep.aida.*;
 import hep.aida.ref.plotter.PlotterRegion;

-import org.lcsim.hps.users.mgraham.alignment.*;

 import hep.physics.matrix.BasicMatrix;
 import hep.physics.matrix.MatrixOp;

-import hep.physics.vec.BasicHep3Matrix;
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Matrix;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;

+import hep.physics.vec.*;

 import java.io.FileWriter;
 import java.io.IOException;
 import java.io.PrintWriter;

-import java.util.*;

+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map;

 import java.util.logging.Level;
 import java.util.logging.Logger;

-import org.lcsim.detector.IDetectorElement;
-import org.lcsim.detector.ITransform3D;
-import org.lcsim.detector.tracker.silicon.ChargeCarrier;

 import org.lcsim.detector.tracker.silicon.SiSensor;

-import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
-import org.lcsim.event.*;

+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;

 import org.lcsim.fit.helicaltrack.*;
 import org.lcsim.hps.event.HPSTransformations;
 import org.lcsim.hps.monitoring.AIDAFrame;
 import org.lcsim.hps.recon.tracking.SvtUtils;
 import org.lcsim.hps.recon.tracking.TrackUtils;
 import org.lcsim.hps.recon.tracking.TrackerHitUtils;

-import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
-import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;

+import org.lcsim.hps.users.mgraham.alignment.RunAlignment;

 import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
 import org.lcsim.recon.tracking.seedtracker.SeedTrack;
 import org.lcsim.util.aida.AIDA;

@@ -56,30 +51,27 @@

 
     private int _nlc = 5;  //the five track parameters
     private int _ngl = 1; //delta(u) and delta(gamma) for each plane

-    //private GlobalParameters _globalParameterList;

     private List<GlobalParameter> _glp = new ArrayList<GlobalParameter>();
     private BasicMatrix _dfdq;

-    //private BasicMatrix _dfdp;

     private HelicalTrackFit _trk;
     private double[] _resid = new double[3];
     private double[] _error = new double[3];
     FileWriter fWriter;
     PrintWriter pWriter;

-    private HPSTransformations _detToTrk;

+    private String _type = "LOCAL";

     boolean _DEBUG=false;
     double smax = 1e3;

+    Hep3Vector _bfield = null;

     ResLimit _resLimits = new ResLimit();

+    private HPSTransformations _detToTrk;

     AlignmentUtils _alignUtils;
     TrackUtils trackUtil;
     TrackerHitUtils trackerHitUtil;

-    Hep3Vector _bfield = null;

     WTrackUtils wutils = new WTrackUtils();
     

-    boolean _doSimHitsResiduals = false;
-    boolean _includeMS = false;

+    private boolean _includeMS = true;

     
     private AIDAFrame plotterFrame;

-    private AIDAFrame plotterFrameSimHit;

     private AIDAFrame plotterFrameSummary;
     private boolean hideFrame = false;
     private AIDA aida = AIDA.defaultInstance();

@@ -99,27 +91,21 @@

-    public MPAlignmentParameters(String outfile, Hep3Vector bfield_vec, boolean debug, boolean hidePlots) {
-        _DEBUG = debug;

+    public MPAlignmentParameters(String outfile) {
+

         _detToTrk = new HPSTransformations();
         _alignUtils = new AlignmentUtils(_DEBUG);
         trackUtil = new TrackUtils();
         trackerHitUtil = new TrackerHitUtils(_DEBUG);

-        hideFrame = hidePlots;
-        _bfield = bfield_vec;
-        this.wutils.setDebug(debug);

+        _bfield = new BasicHep3Vector(0., 0., 1.);
+        this.wutils.setDebug(this._DEBUG);

         try {

-//open things up

             fWriter = new FileWriter(outfile);
             pWriter = new PrintWriter(fWriter);
         } catch (IOException ex) {
             Logger.getLogger(RunAlignment.class.getName()).log(Level.SEVERE, null, ex);
         }
         

-        //_globalParameterList = new GlobalParameters();
-        //if(_DEBUG) _globalParameterList.print();
-        
-        

         makeAlignmentPlots();
         
         

@@ -136,7 +122,29 @@

         _resLimits.add(s,l,d, low,high);
     }
     

-    public void PrintResidualsAndDerivatives(Track track, int itrack, String type, List<SimTrackerHit> simhits) {

+    public void setHideFrame(boolean hide) {
+        this.hideFrame = hide;
+    }
+    
+    public void setDebug(boolean debug) {
+        this._DEBUG = debug;
+    }
+
+    public void setBfield(double bfield) {
+        _bfield = new BasicHep3Vector(0., 0., bfield);
+    }
+    
+    public void setType(String t) {
+        this._type = t;
+    }
+    
+    public void setIncludeMS(boolean include) {
+        this._includeMS = include;
+    }
+
+    
+    
+    public void PrintResidualsAndDerivatives(Track track, int itrack) {

         
         SeedTrack st = (SeedTrack) track;
         SeedCandidate seed = st.getSeedCandidate();

@@ -145,8 +153,8 @@

         List<TrackerHit> hitsOnTrack = track.getTrackerHits();
         String half = hitsOnTrack.get(0).getPosition()[2]>0 ? "top" : "bottom";
         pWriter.printf("TRACK %s (%d)\n",half,itrack);

-        aida.histogram1D("Track Chi2 "+ half).fill(track.getChi2());
-        if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Loop over " + hitsOnTrack.size() + " hitsOnTrack");

+        aida.cloud1D("Track Chi2 "+ half).fill(track.getChi2());
+        if(_DEBUG) System.out.printf("%s: track %d has %d hits\n",this.getClass().getSimpleName(),itrack,hitsOnTrack.size());

         for (TrackerHit hit : hitsOnTrack) {
             
             HelicalTrackHit htc = (HelicalTrackHit) hit;

@@ -164,36 +172,29 @@

 
 
             List<HelicalTrackStrip> clusterlist = cross.getStrips();

-            
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Loop over " + clusterlist.size() + " clusterlist for this hitontrack");
-            
-            for (HelicalTrackStrip cl : clusterlist) {

-                if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": cluster size " + cl.rawhits().size());
-                
-                if(type=="GLOBAL") {
-                    
-                    //CalculateLocalDerivativesGLOBAL(cl);
-                    CalculateGlobalDerivatives(cl);
-                    CalculateResidualGLOBAL(cl, msdrphi, msdz);

+            if(_DEBUG) {
+                System.out.printf("%s: This hit has %d clusters msdrphi=%.4f msdz=%.4f\n",this.getClass().getSimpleName(),clusterlist.size(),msdrphi,msdz);
+            }

-                } else {

+            for (HelicalTrackStrip cl : clusterlist) {
+                if(!"GLOBAL".equals(_type)) {

                     CalculateResidual(cl, msdrphi, msdz);

-                    
-                    if(_doSimHitsResiduals) CalculateResidualSim(cl, msdrphi, msdz,simhits);
-                    

                     CalculateLocalDerivatives(cl);
                     CalculateGlobalDerivatives(cl);                    

-                    PrintStripResidualsNew(cl);

+                    PrintStripResiduals(cl);
+                } else {                    
+                    //CalculateLocalDerivativesGLOBAL(cl);
+                    CalculateGlobalDerivatives(cl);
+                    CalculateResidualGLOBAL(cl, msdrphi, msdz);

                 }
             }
         }
         
         if(itrack%50==0) this.updatePlots();

-        //AddTarget(0.1, 0.02);

     }
 

-    private BasicMatrix CalculateLocalDerivativesOld(HelicalTrackStrip strip) {

+    private BasicMatrix CalculateLocalDerivativesOld(double xint) {

          //get track parameters.
         double d0 = _trk.dca();
         double z0 = _trk.z0();

@@ -201,7 +202,6 @@

         double phi0 = _trk.phi0();
         double R = _trk.R();
 //strip origin is defined in the tracking coordinate system (x=beamline)

-        double xint = strip.origin().x();

         double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
         double phi = -s/R + phi0;
         double[][] dfdq = new double[3][5];

@@ -224,12 +224,10 @@

     }
     
     private void CalculateLocalDerivatives(HelicalTrackStrip strip) {

-        
-        BasicMatrix dfdqGlobalOld = CalculateLocalDerivativesOld(strip);
-        
-        BasicMatrix dfdqGlobal = _alignUtils.calculateLocalHelixDerivatives(_trk, strip, smax, _nlc);
-        Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
-        //_dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal);

+        double xint = strip.origin().x();
+        BasicMatrix dfdqGlobalOld = CalculateLocalDerivativesOld(xint);
+        BasicMatrix dfdqGlobal = _alignUtils.calculateLocalHelixDerivatives(_trk, xint);
+        Hep3Matrix trkToStrip = trackerHitUtil.getTrackToStripRotation(strip);

         _dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal);
         
         if (_DEBUG) {

@@ -240,22 +238,17 @@

             double slope = _trk.slope();
             double phi0 = _trk.phi0();
             double R = _trk.R();

-    //strip origin is defined in the tracking coordinate system (x=beamline)
-            double xint = strip.origin().x();

             double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
             double phi = -s/R + phi0;
             double[] trackpars = {d0, z0, slope, phi0, R, s, xint};

-            System.out.printf("%10s%10s%10s%10s%10s\n","d0","z0","slope","phi0","R");
-            System.out.printf("%10.5f%10.5f%10.5f%10.5f%10.5f\n", d0, z0, slope, phi0, R);
-            System.out.println("Strip Origin: ");
-            System.out.println(strip.origin());
-            System.out.println("s         xint");
-            System.out.printf("%5.5f %5.5f\n", s, xint);
-            System.out.println("Compare local derivatives ");
-            System.out.println("dfdqGlobal:");
-            System.out.println(dfdqGlobal.toString());
-            System.out.println("dfdqGlobalOld:");
-            System.out.println(dfdqGlobalOld.toString());

+            System.out.printf("%s: --- CalculateLocalDerivatives Result --- \n",this.getClass().getSimpleName());
+            System.out.printf("%s: Strip Origin %s \n",this.getClass().getSimpleName(),strip.origin().toString());
+            System.out.printf("%s: %10s%10s%10s%10s%10s%10s%10s\n",this.getClass().getSimpleName(),"d0","z0","slope","phi0","R", "xint", "s");
+            System.out.printf("%s: %10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n",this.getClass().getSimpleName(), d0, z0, slope, phi0, R,xint,s);
+            System.out.printf("%s: Local derivatives:\n",this.getClass().getSimpleName());
+            System.out.printf("%s: \n%s\n",this.getClass().getSimpleName(),dfdqGlobal.toString());
+            System.out.printf("%s: OLD Local derivatives:\n",this.getClass().getSimpleName());
+            System.out.printf("%s: \n%s\n",this.getClass().getSimpleName(),dfdqGlobalOld.toString());

         }
     
     }

@@ -267,30 +260,38 @@

     
     private void CalculateGlobalDerivatives(HelicalTrackStrip strip) {
        

+        //Find interception with plane that the strips belongs to
+        Hep3Vector trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()));
+        
+        if(Double.isNaN(trkpos.x())) {
+            System.out.printf("%s: error this trkpos is wrong %s\n",this.getClass().getSimpleName(),trkpos.toString());
+            System.out.printf("%s: origin %s trk \n%s\n",this.getClass().getSimpleName(),strip.origin(),_trk.toString());
+            trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()),true);
+            System.exit(1);
+        }

         
         double slope = _trk.slope();

-        double xint = strip.origin().x();

+        double xint = trkpos.x();

         double xr = 0.0;
         double yr = 0.0;
         double d0 = _trk.dca();
         double phi0 = _trk.phi0();
         double R = _trk.R();
         double z0 = _trk.z0();

-        double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);

+        double s = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);

         double phi = -s/R + phi0;
         double umeas = strip.umeas();
         Hep3Vector corigin = strip.origin();
         double vmeas = corigin.y(); //THis is wrong

-        

+        int layer = strip.layer();
+        int side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? 10000 : 20000;

 
         if(_DEBUG) {

-            System.out.println("Calculate global derivatives for this strip hit in layer " + strip.layer());
-            System.out.printf(" %10s%10s%10s%10s%10s%10s%10s\n","d0","z0","slope","phi0","R","xint","phi");
-            System.out.printf(" %10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n", d0, z0, slope, phi0, R,xint,phi);
-            

+            System.out.printf("%s: --- Calculate global derivatives ---\n",this.getClass().getSimpleName());
+            System.out.printf("%s: Side %d, layer %d, strip origin %s\n",this.getClass().getSimpleName(),side,layer,corigin.toString());
+            System.out.printf("%s: %10s%10s%10s%10s%10s%10s%10s%10s%10s\n",this.getClass().getSimpleName(),"d0","z0","slope","phi0","R","xint","phi", "xint","s");
+            System.out.printf("%s: %10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n",this.getClass().getSimpleName(), d0, z0, slope, phi0, R,xint,phi,xint,s);

         }

-        if(1==1)
-            return;

         
         //1st index = alignment parameter (only u so far)
         //2nd index = residual coordinate (on du so far)

@@ -316,18 +317,25 @@

         //into the sensor frame because that is where the residuals are calculated
         
        

-        //Which half of the detecotor?
-        int side = strip.origin().z()>0 ? 10000 : 20000;
-        int layer = strip.layer();
-        

+       

        //Rotation matrix from the tracking fram to the sensor/strip frame

-        Hep3Matrix T = this.trackerHitUtil.getTrackToStripRotation(strip);

+        Hep3Matrix T = trackerHitUtil.getTrackToStripRotation(strip);

         
         //Derivatives of the rotation matrix deltaR' w.r.t. rotations a,b,c around axis x,y,z
         Hep3Matrix ddeltaRprime_da = new BasicHep3Matrix(0,0,0,0,0,1,0,-1,0);
         Hep3Matrix ddeltaRprime_db = new BasicHep3Matrix(0,0,-1,0,0,0,1,0,0);
         Hep3Matrix ddeltaRprime_dc = new BasicHep3Matrix(0,1,0,-1,0,0,0,0,0);
 

+        if(_DEBUG) {
+            System.out.printf("%s: Rotation matrx from tracking/global to local frame T:\n %s\n",this.getClass().getSimpleName(),T.toString());
+            System.out.printf("%s: Derivatives of the rotation matrix deltaR' w.r.t. rotation a,b,c around x,y,z axis:\n",this.getClass().getSimpleName());
+            System.out.printf("%s: ddeltaRprime_da:\n %s\n",this.getClass().getSimpleName(),ddeltaRprime_da.toString());
+            System.out.printf("%s: ddeltaRprime_db:\n %s\n",this.getClass().getSimpleName(),ddeltaRprime_db.toString());
+            System.out.printf("%s: ddeltaRprime_dc:\n %s\n",this.getClass().getSimpleName(),ddeltaRprime_dc.toString());
+        }
+        
+        
+

       
         //****************************************************************************
         //Calculate jacobian da/db

@@ -340,14 +348,60 @@

         Hep3Vector dq_dy = VecOp.mult(T, deltaY_gl);
         Hep3Vector dq_dz = VecOp.mult(T, deltaZ_gl);
         

+        if(_DEBUG) {
+            System.out.printf("%s: - Upper left 3x3 of Jacobian da/db: dq_dx,dq_dy,dq_dz\n",this.getClass().getSimpleName());
+            System.out.printf("%s: dq_dx: %s\n",this.getClass().getSimpleName(),dq_dx.toString());
+            System.out.printf("%s: dq_dy: %s\n",this.getClass().getSimpleName(),dq_dy.toString());
+            System.out.printf("%s: dq_dz: %s\n",this.getClass().getSimpleName(),dq_dz.toString());
+        }
+        
+       
+        
+        if(_DEBUG) {
+            System.out.printf("%s: - Upper right 3x3 of Jacobian da/db: dq_dx,dq_dy,dq_dz\n",this.getClass().getSimpleName());
+        }
+        
+        

         // Upper right 3x3

-        Hep3Vector x_vec = new BasicHep3Vector(xint,corigin.y(),corigin.z()); //THIS IS WRONG!!!!
-        Hep3Vector x_vec_tmp = VecOp.mult(ddeltaRprime_da, x_vec);
-        Hep3Vector dq_da = VecOp.mult(T,x_vec_tmp);

+        Hep3Vector x_vec = trkpos; //position around where the small rotation happens 
+        Hep3Vector x_vec_tmp = VecOp.mult(ddeltaRprime_da, x_vec); // derivative of the position
+        Hep3Vector x_vec_tmp2 = VecOp.sub(x_vec_tmp, x_vec); // subtract position
+        Hep3Vector dq_da = VecOp.mult(T,x_vec_tmp2); //rotated into local frame
+        
+        if(_DEBUG) {
+            System.out.printf("%s: position %s rotation derivative w.r.t. a ddeltaR'/da %s\n",this.getClass().getSimpleName(),x_vec.toString(),x_vec_tmp.toString());
+            System.out.printf("%s: subtracted %s and rotated to local %s \n",this.getClass().getSimpleName(),x_vec_tmp2.toString(),dq_da.toString());
+        }
+        
+        

         x_vec_tmp = VecOp.mult(ddeltaRprime_db, x_vec);

-        Hep3Vector dq_db = VecOp.mult(T,x_vec_tmp);

+        x_vec_tmp2 = VecOp.sub(x_vec_tmp, x_vec); 
+        Hep3Vector dq_db = VecOp.mult(T,x_vec_tmp2);
+        
+        if(_DEBUG) {
+            System.out.printf("%s: position %s rotation derivative w.r.t. a ddeltaR'/db %s\n",this.getClass().getSimpleName(),x_vec.toString(),x_vec_tmp.toString());
+            System.out.printf("%s: subtracted %s and rotated to local %s \n",this.getClass().getSimpleName(),x_vec_tmp2.toString(),dq_db.toString());
+        }
+        

         x_vec_tmp = VecOp.mult(ddeltaRprime_dc, x_vec);

-        Hep3Vector dq_dc = VecOp.mult(T,x_vec_tmp);

+        x_vec_tmp2 = VecOp.sub(x_vec_tmp, x_vec); 
+        Hep3Vector dq_dc = VecOp.mult(T,x_vec_tmp2);
+
+        if(_DEBUG) {
+            System.out.printf("%s: position %s rotation derivative w.r.t. a ddeltaR'/dc %s\n",this.getClass().getSimpleName(),x_vec.toString(),x_vec_tmp.toString());
+            System.out.printf("%s: subtracted %s and rotated to local %s \n",this.getClass().getSimpleName(),x_vec_tmp2.toString(),dq_dc.toString());
+        }
+        
+    
+        if(_DEBUG) {
+            System.out.printf("%s: Summary:\n",this.getClass().getSimpleName());
+            System.out.printf("%s: dq_da: %s\n",this.getClass().getSimpleName(),dq_da.toString());
+            System.out.printf("%s: dq_db: %s\n",this.getClass().getSimpleName(),dq_db.toString());
+            System.out.printf("%s: dq_dc: %s\n",this.getClass().getSimpleName(),dq_dc.toString());
+        }
+        
+        
+       

         
         // Lower left 3x3
         //deltaR = (alpha,beta,gamma) the rotation alignment parameters in the local frame

@@ -355,7 +409,19 @@

         Hep3Vector ddeltaR_dy = new BasicHep3Vector(0,0,0);
         Hep3Vector ddeltaR_dz = new BasicHep3Vector(0,0,0);
 

+        
+        if(_DEBUG) {
+            System.out.printf("%s: - Lower left 3x3 of Jacobian ddeltaR/d(x,y,z): dDeltaR_dx,dDeltaR_dy,dDeltaR_dz\n",this.getClass().getSimpleName());
+            System.out.printf("%s: ddeltaR_dx: %s\n",this.getClass().getSimpleName(),ddeltaR_dx.toString());
+            System.out.printf("%s: ddeltaR_dy: %s\n",this.getClass().getSimpleName(),ddeltaR_dy.toString());
+            System.out.printf("%s: ddeltaR_dz: %s\n",this.getClass().getSimpleName(),ddeltaR_dz.toString());
+        }
+        
+        
+     
+        

         // Lower right 3x3

+/*

         //deltaR = (alpha,beta,gamma) the rotation alignment parameters in the local frame
         //Expressing T in Euler angles i,j,k
         double i = 0;

@@ -374,7 +440,15 @@

         Hep3Vector ddeltaR_da = new BasicHep3Vector(dalpha_da,dbeta_da,dgamma_da);
         Hep3Vector ddeltaR_db = new BasicHep3Vector(dalpha_db,dbeta_db,dgamma_db);
         Hep3Vector ddeltaR_dc = new BasicHep3Vector(dalpha_dc,dbeta_dc,dgamma_dc);

+*/
+        
+        if(_DEBUG) {
+            System.out.printf("%s: - Lower right 3x3 of Jacobian ddeltaR/d(a,b,c): \n",this.getClass().getSimpleName());
+            System.out.printf("%s: T: %s\n",this.getClass().getSimpleName(),T.toString());
+        }

+        
+        

         //Now fill the Jacobian 6x6 matrix
         BasicMatrix da_db = new BasicMatrix(6,6);
         //upper left 3x3

@@ -398,15 +472,23 @@

         da_db.setElement(1,5,dq_dc.y());
         da_db.setElement(2,5,dq_dc.z());
         //lower right 3x3

-        da_db.setElement(3,3,ddeltaR_da.x());
-        da_db.setElement(4,3,ddeltaR_da.y());
-        da_db.setElement(5,3,ddeltaR_da.z());
-        da_db.setElement(3,4,ddeltaR_db.x());
-        da_db.setElement(4,4,ddeltaR_db.y());
-        da_db.setElement(5,4,ddeltaR_db.z());
-        da_db.setElement(3,5,ddeltaR_dc.x());
-        da_db.setElement(4,5,ddeltaR_dc.y());
-        da_db.setElement(5,5,ddeltaR_dc.z());

+        da_db.setElement(3,3,T.e(0, 0));
+        da_db.setElement(4,3,T.e(1, 0));
+        da_db.setElement(5,3,T.e(2, 0));
+        da_db.setElement(3,4,T.e(0, 1));
+        da_db.setElement(4,4,T.e(1, 1));
+        da_db.setElement(5,4,T.e(2, 1));
+        da_db.setElement(3,5,T.e(0, 2));
+        da_db.setElement(4,5,T.e(1, 2));
+        da_db.setElement(5,5,T.e(2, 2));
+//        da_db.setElement(4,3,ddeltaR_da.y());
+//        da_db.setElement(5,3,ddeltaR_da.z());
+//        da_db.setElement(3,4,ddeltaR_db.x());
+//        da_db.setElement(4,4,ddeltaR_db.y());
+//        da_db.setElement(5,4,ddeltaR_db.z());
+//        da_db.setElement(3,5,ddeltaR_dc.x());
+//        da_db.setElement(4,5,ddeltaR_dc.y());
+//        da_db.setElement(5,5,ddeltaR_dc.z());

         //lower left 3x3
         da_db.setElement(3,0,ddeltaR_dx.x());
         da_db.setElement(4,0,ddeltaR_dx.y());

@@ -419,18 +501,30 @@

         da_db.setElement(5,2,ddeltaR_dz.z());
         
         if(_DEBUG) {

-            System.out.println("da_db: \n" + da_db.toString());

+            System.out.printf("%s: da_db:\n%s \n",this.getClass().getSimpleName(),da_db.toString());
+            System.out.printf("%s: det(da_db) = %.3f \n",this.getClass().getSimpleName(),da_db.det());

         }
         
         

-        BasicMatrix db_da = da_db;
-        db_da.invert();

+        BasicMatrix db_da = (BasicMatrix) MatrixOp.inverse(da_db);
+        //db_da.invert();

 
         
         if(_DEBUG) {

-            System.out.println("db_da: \n" + db_da.toString());

+            System.out.printf("%s: db_da: \n%s \n",this.getClass().getSimpleName(),db_da.toString());
+            BasicMatrix prod = (BasicMatrix) MatrixOp.mult(db_da,da_db);
+            System.out.printf("%s: db_da*da_db: \n%s \n",this.getClass().getSimpleName(),prod.toString());

         }
 

+        
+        
+        if(1==1)
+            return;
+        
+        
+        
+        
+        

         //****************************************************************************
         // Calcualte the global derivatives in the global frame dz/db
         // where z is the residual in the global frame and b are the alignment 

@@ -664,51 +758,54 @@

     
     
     private void CalculateResidual(HelicalTrackStrip strip, double msdrdphi, double msdz) {

-        if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": CalculateResidual");

+        if(_DEBUG) System.out.printf("%s: --- CalculateResidual ---\n",this.getClass().getSimpleName());

 
         Hep3Vector u = strip.u();
         Hep3Vector v = strip.v();
         Hep3Vector w = strip.w();
         Hep3Vector eta = VecOp.cross(u,v);
         Hep3Vector corigin = strip.origin();

+        String side = SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom";
+
+        double bfield = Math.abs(this._bfield.z());
+        if(_DEBUG) {
+            System.out.printf("%s: Finding interception point for residual calculation (B=%s)\n",this.getClass().getSimpleName(),this._bfield.toString());
+            double xint_wrong = trackUtil.calculateHelixInterceptXPlane(_trk, strip);
+            double s_wrong = HelixUtils.PathToXPlane(_trk, xint_wrong, smax, _nlc).get(0);
+            Hep3Vector trkpos_wrong= HelixUtils.PointOnHelix(_trk, s_wrong);  
+            System.out.printf("%s: using wrong method xint_wrong=%.3f s_wrong=%.3f -> trkpos_wrong=%s\n",this.getClass().getSimpleName(),xint_wrong,s_wrong,trkpos_wrong.toString());
+        }
+        //Find interception with plane that the strips belongs to
+        Hep3Vector trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield);
+        
+        if(_DEBUG) {
+            System.out.printf("%s: found interception point at %s \n",this.getClass().getSimpleName(),trkpos.toString());
+        }
+        
+
+        if(Double.isNaN(trkpos.x()) || Double.isNaN(trkpos.y()) || Double.isNaN(trkpos.z())) {
+            System.out.printf("%s: failed to get interception point (%s) \n",this.getClass().getSimpleName(),trkpos.toString());
+            System.out.printf("%s: track params\n%s\n",this.getClass().getSimpleName(),_trk.toString());
+            System.out.printf("%s: track pT=%.3f chi2=[%.3f][%.3f] \n",this.getClass().getSimpleName(),_trk.pT(bfield),_trk.chisq()[0],_trk.chisq()[1]);
+            trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield,true);
+            System.exit(1);
+        }

-        if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Finding interception point for residual calculation");

-        WTrack wtrack = new WTrack(_trk,Math.abs(_bfield.z()),true); //flip the charge since I'm looking at fitted tracks with B negative?

-        Hep3Vector trkpos = wutils.getHelixAndPlaneIntercept(wtrack, strip.origin(), eta, new BasicHep3Vector(0,0,1));

         double xint = trkpos.x();
         double phi0 = _trk.phi0();
         double R = _trk.R();

-        //double xint = strip.origin().x();

         double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
         double phi = -s/R + phi0;

-        if(_DEBUG) {
-        double xint_wrong = trackUtil.calculateHelixInterceptXPlane(_trk, strip);
-        double s_wrong = HelixUtils.PathToXPlane(_trk, xint_wrong, smax, _nlc).get(0);
-        Hep3Vector trkpos_wrong= HelixUtils.PointOnHelix(_trk, s_wrong);   
-        System.out.println(this.getClass().getSimpleName() + ": trkposdiff " + VecOp.sub(trkpos, trkpos_wrong).toString() + " from trkpos_cor " + trkpos.toString() + " trkpos_wrong " + trkpos_wrong.toString());
-        }

-        //System.out.println("trkpos = "+trkpos.toString());
-        //System.out.println("origin = "+corigin.toString());

+        

         
         Hep3Vector mserr = new BasicHep3Vector(msdrdphi * Math.sin(phi), msdrdphi * Math.sin(phi), msdz);
         Hep3Vector vdiffTrk = VecOp.sub(trkpos, corigin);
         Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
         Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk);

-        //Hep3Vector mserrrot = VecOp.mult(trkToStrip, mserr);
-        //int idiffbin = 0;

-        //double  idiffbin = (vdiffTrk.y())/10;
-//        if(idiffbin<-2) {
-//            idiffbin = -2;
-//            System.out.println("WARNING vdiffTrk.y() = " + vdiffTrk.y() + " merge to -> idiffbin= " + idiffbin);
-//        }
-//        if(idiffbin>2) {
-//            idiffbin = 2;
-//            System.out.println("WARNING vdiffTrk.y() = " + vdiffTrk.y() + " merge to -> idiffbin= " + idiffbin);
-//        }

         
         
         double umc = vdiff.x();

@@ -721,288 +818,38 @@

         double vError = (strip.vmax() - strip.vmin()) / Math.sqrt(12);
         double wmeas = 0;
         double wError = 10.0/Math.sqrt(12); //0.001;

-        //if(idiffbin==0) {
-            _resid[0] = umeas - umc;
-            _error[0] = this._includeMS ? Math.sqrt(uError * uError + msuError * msuError) : uError;
-            _resid[1] = vmeas - vmc;
-            _error[1] = vError;
-            _resid[2] = wmeas - wmc;
-            _error[2] = wError;
-            //if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": idiffbin= " + idiffbin + " resid[0]=" + _resid[0]);
-        //} else {
-        //    _resid[0] = 9999999.9;        
-        //}
-        
-        
-        //Calcualte the distance from the center of the sensor
-        String side = corigin.z()>0. ? "top" : "bottom";
-        //if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": idiffbin= " + idiffbin + " filling ydiff bist with " + _resid[0]);
-        //aida.histogram1D("res_u_vs_ydiff_"+idiffbin+"_layer_" + strip.layer() + "_" + side).fill(_resid[0]);
-        aida.profile1D("res_u_vs_ydiff_layer_" + strip.layer() + "_" + side).fill(vdiffTrk.y(),_resid[0]);
-        
-        if (_DEBUG) {
-            System.out.println("---- " + this.getClass().getSimpleName() + " CalculateResidual ----");
-            System.out.println("Strip Origin: ");
-            System.out.println(corigin.toString());
-            System.out.println("Position on the track (tracking coordinates) at the strip:");
-            System.out.println(trkpos.toString());
-            System.out.println("Difference between the track position and strip origin (tracking coordinates) at the strip:");
-            System.out.println("vdiffTrk :");
-            System.out.println(vdiffTrk.toString());
-            System.out.println("Difference between the track position and strip origin (\"strip\" coordinates) at the strip:");
-            System.out.println("vdiff :");
-            System.out.println(vdiff.toString());
-            System.out.println("u :");
-            System.out.println(u.toString());
-            System.out.println("umeas = " + umeas + ";  umc = " + umc + " ( prediction)");
-            System.out.println("udiff = " + _resid[0] + " +/- " + _error[0] + " ( uError=" + uError + ", msuError=" + msuError + ")");
-            System.out.println("MS: drdphi=" + msdrdphi  + ", dz=" + msdz);
-            System.out.println("MS: phi=" + phi  + " => msvec=" + mserr.toString());
-            System.out.println("MS: msuError = " + msuError + " (msvec*u = " + mserr.toString() + " * " + u.toString() + ")" );
-            //System.out.println("MS: msvec rotated to strip?  = " + mserrrot.toString());
-            
-        }

+        _resid[0] = umeas - umc;
+        _error[0] = this._includeMS ? Math.sqrt(uError * uError + msuError * msuError) : uError;
+        _resid[1] = vmeas - vmc;
+        _error[1] = vError;
+        _resid[2] = wmeas - wmc;
+        _error[2] = wError;

-    }
-    
-    
-     private void CalculateResidualSim(HelicalTrackStrip strip, double msdrdphi, double msdz, List<SimTrackerHit> simTrackerHits) {

-         
-        Hep3Vector u = strip.u();
-        Hep3Vector v = strip.v();
-        Hep3Vector w = strip.w();
-        Hep3Vector corigin = strip.origin();
-        double phi0 = _trk.phi0();
-        double R = _trk.R();
-        //double xint = strip.origin().x();
-        double xint = trackUtil.calculateHelixInterceptXPlane(_trk, strip);
-        double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
-        double phi = -s/R + phi0;
-        Hep3Vector trkpos = HelixUtils.PointOnHelix(_trk, s);

-        //System.out.println("trkpos = "+trkpos.toString());
-        //System.out.println("origin = "+corigin.toString());
-        
-        Hep3Vector mserr = new BasicHep3Vector(msdrdphi * Math.sin(phi), msdrdphi * Math.sin(phi), msdz);
-        Hep3Vector vdiffTrk = VecOp.sub(trkpos, corigin);
-        Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
-        Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk);
-        //Hep3Vector mserrrot = VecOp.mult(trkToStrip, mserr);
-        double umc = vdiff.x();
-        double vmc = vdiff.y();
-        double wmc = vdiff.z();
-        double umeas = strip.umeas();
-        double uError = strip.du();
-        double msuError = VecOp.dot(mserr, u);
-        double vmeas = 0;
-        double vError = (strip.vmax() - strip.vmin()) / Math.sqrt(12);
-        double wmeas = 0;
-        double wError = 10.0/Math.sqrt(12); //0.001;
-//        //System.out.println("strip error="+uError+"; ms error ="+msuError);
-//        _resid[0] = umeas - umc;
-//        _error[0] = Math.sqrt(uError * uError + msuError * msuError);
-//        _resid[1] = vmeas - vmc;
-//        _error[1] = vError;
-//        _resid[2] = wmeas - wmc;
-//        _error[2] = wError;
-        if (_DEBUG) {
-            System.out.println("---- " + this.getClass().getSimpleName() + " CalculateResidualSim ----");
-            System.out.println("Strip Origin: ");
-            System.out.println(corigin.toString());
-            System.out.println("Position on the track (tracking coordinates) at the strip:");
-            System.out.println(trkpos.toString());
-            System.out.println("Difference between the track position and strip origin (tracking coordinates) at the strip:");
-            System.out.println("vdiffTrk :");
-            System.out.println(vdiffTrk.toString());
-            System.out.println("Difference between the track position and strip origin (\"strip\" coordinates) at the strip:");
-            System.out.println("vdiff :");
-            System.out.println(vdiff.toString());
-            System.out.println("u :");
-            System.out.println(u.toString());
-            System.out.println("umeas = " + umeas + ";  umc = " + umc + " ( prediction)");
-            System.out.println("udiff = " + _resid[0] + " +/- " + _error[0] + " ( uError=" + uError + ", msuError=" + msuError + ")");
-            System.out.println("MS: drdphi=" + msdrdphi  + ", dz=" + msdz);
-            System.out.println("MS: phi=" + phi  + " => msvec=" + mserr.toString());
-            System.out.println("MS: msuError = " + msuError + " (msvec*u = " + mserr.toString() + " * " + u.toString() + ")" );
-            
-        }
-        
-        if(simTrackerHits!=null && simTrackerHits.size()>0) {
-            List<SimTrackerHit> simHits = trackerHitUtil.stripClusterToSimHits(strip, simTrackerHits,true);
-            if(simHits.size()>0) {
-                //Should only be one!
-                if(simHits.size()>1) {
-                    System.out.println(this.getClass().getSimpleName() + " ERROR multiple simhits for this strip");
-                    System.exit(1);
-                }
-                if(_DEBUG) System.out.println( simHits.size() + " sim hits for this raw hit at: ");
-                SimTrackerHit simHit =simHits.get(0);
-                if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in JLab frame" + simHit.getPositionVec().toString());
-                //Transformation between the JLAB and tracking coordinate systems
-                Hep3Matrix detToTrackMatrix = trackerHitUtil.detToTrackRotationMatrix();
-                Hep3Vector simHitPos = VecOp.mult(detToTrackMatrix, simHit.getPositionVec());
-                if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in tracking frame" + simHitPos.toString());
-
-                String side = strip.origin().z()>0 ? "top" : "bottom";
-                Hep3Vector stripPosition = trackerHitUtil.getClusterPosition(strip,true);
-                aida.histogram1D("simhitstripres_z_layer" + strip.layer() + "_" + side).fill(simHitPos.z()-stripPosition.z());
-                aida.histogram1D("simhitstripres_y_layer" + strip.layer() + "_" + side).fill(simHitPos.y()-stripPosition.y());
-                aida.histogram1D("simhitstripres_x_layer" + strip.layer() + "_" + side).fill(simHitPos.x()-stripPosition.x());
-
-                //Transform the sim hit to do a u residual
-                Hep3Vector vDiffSimTrk = VecOp.sub(simHitPos, corigin);
-                //Rotate this into the local sensor frame
-                Hep3Vector vdiffSim = VecOp.mult(trkToStrip, vDiffSimTrk);
-                double umc_sim = vdiffSim.x();
-                double vmc_sim = vdiffSim.y();
-                double wmc_sim = vdiffSim.z();
-                aida.histogram1D("simhitstripres_u_layer" + strip.layer() + "_" + side).fill(umeas - umc_sim);
-                aida.histogram1D("simhitstripres_v_layer" + strip.layer() + "_" + side).fill(vmeas - vmc_sim);
-                aida.histogram1D("simhitstripres_w_layer" + strip.layer() + "_" + side).fill(wmeas - wmc_sim);
-            
-            }
-            
-        } else {
-            if(_DEBUG) System.out.println("No simTrackerHit collection");
-        }
-           
-            
-        

+        //Fill residuals vs distrance from center of plane in the v directions
+        aida.profile1D("res_u_vs_ydiff_layer_" + strip.layer() + "_" + side).fill(vdiffTrk.y(),_resid[0]);

-    }
-    
-     
-     
-     public void CalculateResidualSimHits(List<HelicalTrackStrip> strips, List<SiTrackerHit> strip1dhits,  List<SimTrackerHit> simTrackerHits) {
-//        SeedTrack st = (SeedTrack) track;
-//        SeedCandidate seed = st.getSeedCandidate();
-//        Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap();
-//        HelicalTrackFit trk = seed.getHelix();
-//         double msdrdphi = 0;
-         double msdz = 0;
-         
-        if (_DEBUG) {
-          System.out.println("---- " + this.getClass().getSimpleName() + " CalculateResidualSimHits ----");
-        }
-        if(strips==null) {
-            if (_DEBUG) {
-                System.out.println(this.getClass().getSimpleName() + ": list of strips is null");
-            }
-            return;
-        }
-         

         if (_DEBUG) {

-          System.out.println(this.getClass().getSimpleName() + ": " + strips.size() + " strips");

+            System.out.printf("%s: ----  CalculateResidual ----\n",this.getClass().getSimpleName());
+            System.out.printf("%s: Strip Origin: %s\n",this.getClass(),corigin.toString());
+            System.out.printf("%s: Position on the track (tracking coordinates) at the strip: %s\n",this.getClass(),trkpos.toString());
+            System.out.printf("%s: vdiffTrk %s\n",this.getClass(),vdiffTrk.toString());
+            System.out.printf("%s: vdiff %s\n",this.getClass(),vdiff.toString());
+            System.out.printf("%s: u %s\n",this.getClass(),u.toString());
+            System.out.printf("%s: umeas = %.4f  umc = %.4f\n",this.getClass(),umeas,umc);
+            System.out.printf("%s: udiff = %.4f +/- %.4f  ( uError=%.4f , msuError=%.4f\n",this.getClass(),_resid[0],_error[0],uError,msuError);
+            System.out.printf("%s: MS: drdphi=%.4f, dz=%.4f\n",this.getClass(),msdrdphi,msdz);
+            System.out.printf("%s: MS: phi=%.4f => msvec=%s\n",this.getClass(),phi,mserr.toString());
+            System.out.printf("%s: MS: msuError = %.4f (msvec*u = %s * %s\n",this.getClass(),msuError,mserr.toString(),u.toString());          

         }

-         for(HelicalTrackStrip strip : strips) {
-             //Rotate into tracking frame
-             Hep3Matrix detToTrackMatrix = this.trackerHitUtil.detToTrackRotationMatrix();
-            
-            Hep3Vector u = VecOp.mult(detToTrackMatrix,strip.u());
-            Hep3Vector v = VecOp.mult(detToTrackMatrix,strip.v());
-            Hep3Vector w = VecOp.mult(detToTrackMatrix,strip.w());
-            Hep3Vector corigin = VecOp.mult(detToTrackMatrix,strip.origin());
-            double umeas = strip.umeas();
-            double vmeas = 0;
-            double wmeas = 0;
-            if (_DEBUG) {
-                System.out.println("----");
-                System.out.println("Strip Origin: ");
-                System.out.println(corigin.toString());
-                System.out.println("SIMHITS");
-            }
-            if(simTrackerHits!=null && simTrackerHits.size()>0) {
-                List<SimTrackerHit> simHits = this.trackerHitUtil.stripClusterToSimHits(strip, simTrackerHits,false);
-                if(simHits.size()>0) {
-                    //Should only be one!
-                    if(simHits.size()>1) {
-                        System.out.println(this.getClass().getSimpleName() + " ERROR multiple simhits for this strip");
-                        System.exit(1);
-                    }
-                       
-                    if(_DEBUG) System.out.println( simHits.size() + " sim hits for this raw hit at: ");
-                    SimTrackerHit simHit =simHits.get(0);
-                    if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in JLab frame" + simHit.getPositionVec().toString());
-                    //Transformation between the JLAB and tracking coordinate systems
-                    Hep3Vector simHitPos = VecOp.mult(detToTrackMatrix, simHit.getPositionVec());
-                    if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in tracking frame" + simHitPos.toString());
-
-                    String side = corigin.z()>0 ? "top" : "bottom";
-                    Hep3Vector stripPosition = this.trackerHitUtil.getClusterPosition(strip,false);
-                    aida.histogram1D("simhitallstripres_z_layer" + strip.layer() + "_" + side).fill(simHitPos.z()-stripPosition.z());
-                    aida.histogram1D("simhitallstripres_y_layer" + strip.layer() + "_" + side).fill(simHitPos.y()-stripPosition.y());
-                    aida.histogram1D("simhitallstripres_x_layer" + strip.layer() + "_" + side).fill(simHitPos.x()-stripPosition.x());
-
-    //                //Transform the sim hit to do a u residual
-    //                Hep3Vector vDiffSimTrk = VecOp.sub(simHitPos, corigin);
-    //                //Rotate this into the local sensor frame
-    //                Hep3Vector vdiffSim = VecOp.mult(trkToStrip, vDiffSimTrk);
-    //                double umc_sim = vdiffSim.x();
-    //                double vmc_sim = vdiffSim.y();
-    //                double wmc_sim = vdiffSim.z();
-    //                aida.histogram1D("simallhitstripres_u_layer" + strip.layer() + "_" + side).fill(umeas - umc_sim);
-    //                aida.histogram1D("simallhitstripres_v_layer" + strip.layer() + "_" + side).fill(vmeas - vmc_sim);
-    //                aida.histogram1D("simallhitstripres_w_layer" + strip.layer() + "_" + side).fill(wmeas - wmc_sim);
-                } 
-
-
-            } else {
-                if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " No simTrackerHit collection");
-            }
-            
-            /*
-            if(strip1dhits!=null && strip1dhits.size()>0) {
-                List<SiTrackerHit> siHits = this.trackerHitUtil.stripClusterToSiHits(strip, strip1dhits, false);
-                if(siHits.size()>0) {
-                    //Should only be one!
-                    if(siHits.size()>1) {
-                        System.out.println(this.getClass().getSimpleName() + " ERROR multiple simhits for this strip");
-                        System.exit(1);
-                    }
-                    SiTrackerHitStrip1D h = (SiTrackerHitStrip1D) siHits.get(0);
-                    if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " SiTrackerHitStrip1D position in JLab frame" + h.getPositionAsVector().toString());
-                    //Transformation between the JLAB and tracking coordinate systems
-                    Hep3Vector siHitPos = VecOp.mult(detToTrackMatrix, h.getPositionAsVector());
-                    if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in tracking frame" + siHitPos.toString());
-
-//                    String side = corigin.z()>0 ? "top" : "bottom";
-//                    Hep3Vector stripPosition = this.getClusterPosition(strip,false);
-//                    aida.histogram1D("simhitallstripres_z_layer" + strip.layer() + "_" + side).fill(simHitPos.z()-stripPosition.z());
-//                    aida.histogram1D("simhitallstripres_y_layer" + strip.layer() + "_" + side).fill(simHitPos.y()-stripPosition.y());
-//                    aida.histogram1D("simhitallstripres_x_layer" + strip.layer() + "_" + side).fill(simHitPos.x()-stripPosition.x());
-
-    //                //Transform the sim hit to do a u residual
-    //                Hep3Vector vDiffSimTrk = VecOp.sub(simHitPos, corigin);
-    //                //Rotate this into the local sensor frame
-    //                Hep3Vector vdiffSim = VecOp.mult(trkToStrip, vDiffSimTrk);
-    //                double umc_sim = vdiffSim.x();
-    //                double vmc_sim = vdiffSim.y();
-    //                double wmc_sim = vdiffSim.z();
-    //                aida.histogram1D("simallhitstripres_u_layer" + strip.layer() + "_" + side).fill(umeas - umc_sim);
-    //                aida.histogram1D("simallhitstripres_v_layer" + strip.layer() + "_" + side).fill(vmeas - vmc_sim);
-    //                aida.histogram1D("simallhitstripres_w_layer" + strip.layer() + "_" + side).fill(wmeas - wmc_sim);
-                } else {
-                        System.out.println(this.getClass().getSimpleName() + " ERROR no SiTrackerHitStrip1D for this strip");
-                        //System.exit(1);
-
-                }
-
-
-            } else {
-                if(_DEBUG) System.out.println("No SiTrackerHitStrip1D matched to this strip");
-            }
-            
-            */
-            
-         }
-            
-        

 
     }
     
     
     

+     
+     

@@ -1042,20 +889,12 @@

-    private void PrintStripResidualsNew(HelicalTrackStrip strip) {

+    private void PrintStripResiduals(HelicalTrackStrip strip) {

         if (_DEBUG) {

-            System.out.println(this.getClass().getSimpleName() + ": PrintStripResidualsNew");
-            int s = 1;
-            if(strip.origin().z()>0) s = 0;
-            System.out.printf("Strip Layer =  %4d\n", strip.layer());
-            System.out.printf("Residuals (u,v,w) : %5.5e %5.5e %5.5e   (limits: [%.1e,%.1e] [%.1e,%.1e] [%.1e,%.1e])\n"
-                    , _resid[0], _resid[1], _resid[2],
-                    _resLimits.getMin(s,strip.layer(),0),
-                    _resLimits.getMax(s,strip.layer(),0),
-                    _resLimits.getMin(s,strip.layer(),1),
-                    _resLimits.getMax(s,strip.layer(),1),
-                    _resLimits.getMin(s,strip.layer(),2),
-                    _resLimits.getMax(s,strip.layer(),2));

+            System.out.println(this.getClass().getSimpleName() + ": --- PrintStripResiduals ---");
+            String side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom";
+            System.out.printf("%s: Strip layer %4d in %s  at origin %s\n",this.getClass().getSimpleName(), strip.layer(), side, strip.origin().toString());
+            System.out.printf("Residuals (u,v,w) : %5.5e %5.5e %5.5e\n", _resid[0], _resid[1], _resid[2]);

             System.out.printf("Errors (u,v,w)    : %5.5e %5.5e %5.5e\n", _error[0], _error[1], _error[2]);
             String[] q = {"d0", "z0", "slope", "phi0", "R"};
             System.out.println("track parameter derivatives");

@@ -1068,7 +907,7 @@

                 System.out.printf("%s  %5.5e %5.5e %5.5e   %5d  %10s\n", "", gl.dfdp(0), gl.dfdp(1), gl.dfdp(2), gl.getLabel(),gl.getName());
                 //System.out.printf("%s  %5.5e %5.5e %5.5e   %5d\n", p[j], _dfdp.e(0, j), _dfdp.e(1, j), _dfdp.e(2, j), _globalLabel[j]);
             }

-

+            System.out.println(this.getClass().getSimpleName() + ": --- END PrintStripResiduals ---");

         }
         pWriter.printf("%d\n", strip.layer());
         // Loop over the three directions u,v,w and print residuals and derivatives

@@ -1257,50 +1096,46 @@

-         aida.tree().cd("/");

+        aida.tree().cd("/");

         plotterFrame = new AIDAFrame();
         plotterFrame.setTitle("Residuals");

-

         plotterFrameSummary = new AIDAFrame();
         plotterFrameSummary.setTitle("Summary");
               

-        plotterFrameSimHit = new AIDAFrame();
-        plotterFrameSimHit.setTitle("Residuals SimHits");
-

         String[] side = {"top","bottom"};
 
         IPlotter plotter_chi2 = af.createPlotterFactory().create();
         plotter_chi2.createRegions(2,1,0);
         plotter_chi2.setTitle("Track Chi2");
         plotterFrame.addPlotter(plotter_chi2);

-        IHistogram hchi2_top = aida.histogram1D("Track Chi2 top",50,0,15);
-        IHistogram hchi2_bot = aida.histogram1D("Track Chi2 bottom",50,0,3);

+        ICloud1D hchi2_top = aida.cloud1D("Track Chi2 top");
+        ICloud1D hchi2_bot = aida.cloud1D("Track Chi2 bottom");

         plotter_chi2.region(0).plot(hchi2_top);
         plotter_chi2.region(1).plot(hchi2_bot);
         
         
         String[] direction = {"u"};
         double xbins_u_res[][] = {

-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01},
-            {-0.01,0.01}
-//                {-0.3,0.3},
-//                {-0.7,0.7},
-//                {-0.7,0.7},
-//                {-0.5,0.5},
-//                {-0.5,0.5},
-//                {-0.5,0.5},
-//                {-2,2},
-//                {-2,2},
-//                {-2,2},
-//                {-2,2}            

+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01},
+//            {-0.01,0.01}
+                {-0.15,0.15},
+                {-0.4,0.4},
+                {-0.7,0.7},
+                {-0.5,0.5},
+                {-0.5,0.5},
+                {-0.5,0.5},

[truncated at 1000 lines; 150 more skipped]