lcd/hep/lcd/recon/tracking/ccd

diff -N CCDECSpecifications.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ CCDECSpecifications.java	10 Aug 2005 20:28:16 -0000	1.1
@@ -0,0 +1,122 @@

+/* CCDSpecifications.java
+
+		Specifies CCD dimensions and performance.
+
+ Created by  Mike Ronan  Jan '99   Used in passing tracker specifications to TrkCandidate
+                                   in recalculating track parameters.
+ Modified by
+
+ */
+
+package hep.lcd.recon.tracking.ccd;
+
+import hep.lcd.recon.tracking.*;
+
+import hep.lcd.event.*;
+import java.io.Serializable;
+
+import hep.lcd.geometry.Detector;
+import hep.lcd.geometry.PropertySet;
+import hep.lcd.geometry.component.*;
+
+final public class CCDECSpecifications extends BaseSpecifications
+{
+    final int MaxLayers = 10;
+	boolean haveEC = false;
+
+    public CCDECSpecifications()
+    {
+     this(false);
+	 System.out.println("Creating CCD EC specifications");	
+    }
+
+    public CCDECSpecifications(boolean deb)
+    {
+      Type = "CCDEC";
+      this.setSpecs(deb);
+    }
+	
+	public boolean haveEndCap()
+	{
+	 return haveEC;	
+	}
+    
+    public void setSpecs(boolean deb)
+    {
+      debug = deb;
+      findType();
+      LayRadii = new double[MaxLayers];
+      LrsHalfLength = new double[MaxLayers];
+      LayOutRad = new double[MaxLayers];
+      for(int i=0; i<MaxLayers; i++)
+      {
+       LayRadii[i] = 0.;
+       LrsHalfLength[i] = 0.;
+       LayOutRad[i] = 0.;
+      }
+      NLayers = 0;
+	  System.out.println("CCD EC specs: Detector is: "+d);
+      PropertySet vpr = (PropertySet) d.get("MultiLayer_VERTX_EndCap");
+	  if(vpr==null) System.out.println("No Vertex Endcap !");
+      MaterialList mlst = (MaterialList) d.get("VERTX_EndCap_Material");
+	  if(mlst == null) System.out.println("No material list for VXD EC");
+      double actof = 0.;
+	  if(mlst != null)
+	  {
+		System.out.println("VXD EC: Processing material list ");  
+       for(int i=0; i<mlst.getNMaterials(); i++)
+       {
+        Material mat = mlst.getMaterial(i);
+        if(!mat.isActive()) actof+=mat.getThickness();
+        else break;
+       }
+      }   
+      double noe = 0.;
+      if(vpr != null)
+	  {
+       noe = vpr.getDouble("num_of_endcaps");
+	   haveEC = true;
+      }
+	  System.out.println("VXD EC has "+noe+" endcap lrs");
+      if(noe > 0.5)
+      { 
+       InnerRadius = vpr.getDouble("inner_radius1");
+       double nl1 = vpr.getDouble("layers1");
+       int inob = (int) noe;
+       int inl1 = (int) nl1;
+       double lthick = 0.03;
+       double zstrt = 0.;
+       double zend = 0.;
+       for(int i =1; i<inob+1; i++)
+       {
+        double strad = vpr.getDouble("inner_radius"+i);
+        double endrad = vpr.getDouble("outer_radius"+i);
+        double nlpb = vpr.getDouble("layers"+i);
+        double zin = vpr.getDouble("z_inner"+i);
+        double zout = vpr.getDouble("z_outer"+i);
+        if(i==1) zstrt = zin;
+        if(i==inob) zend = zout;
+        lthick = zout-zin;
+        lthick/=nlpb;
+        int inlpb = (int) nlpb;
+        for(int j=0; j<inlpb; j++)
+        {
+         if(NLayers > MaxLayers-1)
+         {
+          System.out.println("CCDECSpecifications::error in number of layers in CCD EC!");
+          NLayers=MaxLayers-1;
+         }  
+         LayRadii[NLayers]=strad;
+         LayOutRad[NLayers] = endrad;
+         LrsHalfLength[NLayers] = zin + actof + j*lthick;
+         NLayers++;
+        }
+       }
+       HalfLength = zend - zstrt;
+       if(debug) System.out.println("Specifications: CCD EC noec="+inob+" total layers="+NLayers); 
+       OuterRadius = vpr.getDouble("outer_radius"+inob);
+      } 	  
+     debug = false;
+    }
+}
+

lcd/hep/lcd/recon/tracking/ccd

diff -u -r1.2 -r1.3
--- VXDFindStrategy.java	20 Jul 2005 17:11:35 -0000	1.2
+++ VXDFindStrategy.java	10 Aug 2005 20:28:16 -0000	1.3
@@ -22,11 +22,16 @@

  public VXDFindStrategy()
  {
      int[] Levels = {2,1};

-     int[] Patterns = {3,3};
-     double[][] Radius = {{40.,15.},{6.,6.}};
-     int[][] pnts = {{90,20},{6,6}};
-     int[][][] tri = {{{97,45,4},{25,11,3},{10,5,2}},{{4,2,0},{4,3,1},{3,2,0}}};

+     int[] Patterns = {2,13};
+     double[][] Radius = {{40.,15.},{8.,8.}};
+     int[][] pnts = {{6,6},{6,6}};
+     int trip[][][] = {{{4,2,0},{4,3,1},{5,2,1},{7,6,5},{8,7,5},{0,0,0},{0,0,0},
+		 {0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0}},
+	     {{3,1,0},{4,2,0},{4,3,1},{5,2,0},{5,3,1},{5,1,0},{6,2,0},{6,5,1},
+	     {6,1,0},{7,5,0},{7,6,1},{8,6,0},{8,7,5}}};
+	 rTol[0] = 0.5;

 	 rTol[1] = 0.5;

+	 zTol[0] = 0.5;

 	 zTol[1] = 0.5;
      for(int i=0; i<2; i++)
      {

@@ -41,10 +46,35 @@

       {
        for(int n=0; n<3; n++)
        {

-        triplets[i][j][n] = tri[i][j][n];

+        triplets[i][j][n] = trip[i][j][n];

        }
       }
      }

- }
-
-}

\ No newline at end of file

+   }
+   public void setNpntsL1(int[] np)
+   {
+	   for(int i=0; i<2; i++)
+	   {
+		 nPtsReq[i][0] = np[i];  
+	   }   
+   }   
+   public void setNpntsL2(int[] np)
+   {
+	   for(int i=0; i<2; i++)
+	   {
+		 nPtsReq[i][1] = np[i];  
+	   }   
+   }   
+   public void setRTol(double[] rt)
+   {
+	 rTol[0]=rt[0];
+	 rTol[1]=rt[1];
+	 return;
+   }
+   public void setZTol(double[] zt)
+   {
+	 zTol[0]=zt[0];
+	 zTol[1]=zt[1];
+   }   
+   
+}

lcd/hep/lcd/recon/tracking/ccd

diff -u -r1.4 -r1.5
--- VXDPatFind.java	20 Jul 2005 17:11:35 -0000	1.4
+++ VXDPatFind.java	10 Aug 2005 20:28:16 -0000	1.5
@@ -13,7 +13,6 @@

 	There is no flag to tell VXDPatFinder whether to use Endcap data or not -
 	it uses whatever hits are available. But VXDReco may remove Endcap hits
 	from hit data arrays passed to VXDPatFinder if Endcap should not be used.

-	Usage of VXD disks is not implemented yet in this version of pattern finder.

 
  Created by Nick Sinev, June 2004
 

@@ -22,6 +21,7 @@

 package hep.lcd.recon.tracking.ccd;
 
 import hep.lcd.recon.tracking.ccd.CCDSpecifications;

+import hep.lcd.recon.tracking.ccd.CCDECSpecifications;

 import hep.lcd.recon.tracking.tpc.TPCSpecifications;
 import hep.lcd.recon.tracking.*;
 import hep.lcd.recon.tracking.ccd.CCDReco;

@@ -47,15 +47,19 @@

     Specifications vspecs;
     Specifications specs;
     Specifications ecspecs;

+	CCDECSpecifications vespecs;

     ThreePntCircle tpcir; 
     CCDReco vxd;
 	//.................................................................
 	// these are parameters to trace reconstruction of particular
 	// track for debugging purposes

+	private boolean traceAll = false;

     private boolean trace = false;
     private int tracehc = 0;

+	boolean printall = true;

     private boolean dotrace = false;
 	private double[] ctrace = {0.,0.,0.};

+	int ntracepr = 0;

 	MCParticle mcptot;
 	MCParticle mcpfcth;
 	MCParticle mclstvx;

@@ -66,12 +70,16 @@

 	//.........................................................................
 	//   parameters set in the setParameters function, extracted from detector
 	//   specifications and constructor inputs

+	boolean useEC = false;
+	boolean useVXEC = false;

     int maxpnts = 0;  // maximum number of hits in all trackers (VXD,CT,EC) to handle
 	boolean dohist = false; //do histogramming ?
     double tdInnerR = 0.;  // central tracker inner radius
     double tdOuterR = 0.;  // central tracker outer radius
     double ecresol = 0.;   // endcap RPhi resolution
     int nvxlrs = 0;        // number of layers in Vertex detector

+	int nvxblrs = 0;       // number of layers in barrel Vertex detector
+	int nvxelrs = 0;       // number of layers in EC Vertex detector

     int nctlrs = 0;        // number of layers in Central Tracker
     int neclrs = 0;        // number of layers in Endcap tracker
     int nltot=20;          // total number of layers

@@ -84,19 +92,21 @@

 	double vxterr = 0.0; //Vertex detector spatial resolution
 	// these are not defined in any ini file. in principle they should.
     double ecrres = 10.;    // the length of endcap radial strip segment 

-    double ecper = 0.3;
-    double ecrer = 0.3;     //search distanse for EC hit

+    double ecper = 0.5;
+    double ecrer = 0.5;     //search distanse for EC hit

     int ecnsegs = 16;       // number of endcap strip segments
 	//..............................................................................
 	// these two values are coefficients for calculating multipple scatterring angle
 	// from curvature radius for SiD. In principle, they should be calculated
 	// from detector specification, but hardcoded here for now. Theta0 = coeff/rcurv

-	double vxdmsco = 0.0228;

+//	double vxdmsco = 0.0228;
+	double vxdmsco = 0.0114;

 	double trkmsco = 0.1437;
 	//..............................................................................
 //    double ctresol = 0.001;
 //    double xyerv = 0.002;
 // These are working variables, not constants

+    int dtyp = 0;

 	double mindistance = 0.;
 	double fcthndis = 0.;
 	double fcthedev = 0.;

@@ -114,7 +124,10 @@

     double ecrersq = 0.;
 	double xclhit = 0.;
 	double yclhit = 0.;

+	double[] vdr;
+    int nexthts = 0;

 	int nancor = 0;

+	int seclhin = -1; 

 //    double msctrm = 0.001;
     double Rcu = 0.;
 //    double minRad = 8.;

@@ -135,26 +148,52 @@

     double[] p2 = new double[2];
     double[] p3 = new double[2];
 	double[] p4 = new double[2];

+    int[] vxdhin = null;

 	int ndegfrvtx = 0;
 	int lstbrlrr = 0;
 	double asshchsq = 0.;
 	int nastrhits = 0;

+	int lstvxhin = 0;
+	int outerLayer  = 0;
+	int middleLayer = 0;
+	int innerLayer  = 0;
+	boolean fromIP = false;

 	TrackerHit[] Thits;
 //	Histogram2D h101 = null;

-//	Histogram h010 = null;
-//	Histogram h011 = null;
-//	Histogram h012 = null;
-//	Histogram h013 = null;
-//	Histogram h014 = null;
-//	Histogram h015 = null;
-//	Histogram h016 = null;
-//	Histogram h017 = null;
-//	Histogram h018 = null;
-//	Histogram h019 = null;
-//	Histogram h021 = null;
-//	Histogram h022 = null;
-//	Histogram h023 = null;
-//	Histogram h024 = null;

+/*	Histogram h010 = null;
+	Histogram h4 = null;
+	Histogram h010_0 = null;
+	Histogram h010_1 = null;
+	Histogram h010_2 = null;
+	Histogram h010_3 = null;
+	Histogram h010_4 = null;
+	Histogram h010_5 = null;
+	Histogram h010_6 = null;
+	Histogram h010_7 = null;
+	Histogram h010_8 = null;
+	Histogram h011 = null;
+	Histogram h012 = null;
+	Histogram h012a = null;
+	Histogram h012d = null;
+	Histogram h012am = null;
+	Histogram h012az = null;
+	Histogram h012nc = null;
+	Histogram2D h112 = null;
+	Histogram h012b = null;
+	Histogram h013 = null;
+	Histogram h014 = null;
+	Histogram h015 = null;
+	Histogram h016 = null;
+	Histogram h017 = null;
+	Histogram h018 = null;
+	Histogram h019 = null;
+	Histogram h021 = null;
+	Histogram h022 = null;
+	Histogram h023 = null;
+	Histogram h024 = null;
+	Histogram h036 = null;
+	Histogram h037 = null;
+	Histogram h038 = null; */

 //	Histogram2D h102 = null;
     DecimalFormat df = new DecimalFormat();
 	// Some constants

@@ -162,6 +201,7 @@

 	double halfPI = 0.5 * Math.PI;
 	double quoterPI = 0.25 * Math.PI;
 	double eighthPI = 0.125 *Math.PI;

+	boolean sameMCPhts = false;

 
     /**  Construct VXD TrackFinder (MaxTracks=400,MaxTrkPts=30). */
     public VXDPatFind()

@@ -174,7 +214,7 @@

     }
     public VXDPatFind(boolean hist,VXDFindStrategy strat)
     {

-   	super(400,30,hist); 

+   	super(400,160,hist); 

         Finder = "VXDPatFind";
         tpcir = new ThreePntCircle();
         df.setMaximumFractionDigits(5);

@@ -186,12 +226,20 @@

         this.Strategy = strat;
     }
 	

+	

 	public void setTrace(int hind)
 	{
 		tracehc = hind;
 		trace = true;
 	}
 	

+	public void setTraceAll(boolean t)
+	{
+		traceAll = t;
+		trace = t;
+		debug = t;
+	}
+	

 	public void setTrace(double[] tc)
 	{
 	  tracehc = -1;

@@ -215,7 +263,7 @@

        }
        if(this.getTracker() == null)
        {

-        System.out.println("TRACKER IS NOT set!");

+        System.err.println("TRACKER IS NOT set!");

         tracker = null;
        } 
        if(this.getTracker() != null)

@@ -223,16 +271,27 @@

         specs = (Specifications) this.getTracker().getSpecifications();
         ecspecs = (Specifications) ((VXDReco)this.getTracker()).getECSpecifications();
         vspecs = (Specifications) ((VXDReco)this.getTracker()).getVxdSpecifications();

+		vespecs = (CCDECSpecifications) ((VXDReco)this.getTracker()).getVxdECSpecifications();

         maxpnts = ((VXDReco)this.getTracker()).getMaxPoints();
 		MaxPoints = maxpnts;

-		System.out.println("Tracker can process "+maxpnts+" hits ");

+//		System.out.println("Tracker can process "+maxpnts+" hits ");
+		useEC=((VXDReco)tracker).includeEC;
+		useVXEC = ((VXDReco)tracker).includeVXDEC;
+//		if(useVXEC) System.out.println("Using Endcap VXD !");

 		mainco = new double[maxpnts];

-        int dtyp=specs.getDetectorType();

+        dtyp=specs.getDetectorType();

         tdInnerR = specs.getInnerRadius();
         tdOuterR = specs.getOuterRadius();
        // ctresol = specs.getRPhiResolution();
         ecresol = ecspecs.getRPhiResolution();

-        nvxlrs = vspecs.getNLayers();

+        nvxblrs = vspecs.getNLayers();
+		vdr = new double[nvxblrs];
+		for(int i=0; i<nvxblrs; i++)
+		 vdr[i]=vspecs.getLayerRadius(i);
+		nvxelrs = vespecs.getNLayers();
+		nvxlrs = nvxblrs;
+		if(useVXEC) nvxlrs = nvxblrs+nvxelrs;
+        vxdhin = new int[nvxlrs];

         nctlrs = specs.getNLayers();
         neclrs = ecspecs.getNLayers();
         nltot = nvxlrs+nctlrs+neclrs;

@@ -241,28 +300,36 @@

         ispin = new int[maxpnts]; 
         ctrInnerR = specs.getLayerRadius(0);
         ctrInnerZ = specs.getLayerZ(0);          

-        int MaxLevel = Strategy.getLevels(dtyp) - 1;
-        double[] minRadius = Strategy.getMinRadius(dtyp);
-        int Npatterns = Strategy.getNPatterns(dtyp);
-        int[][] triplets = Strategy.getTriplets(dtyp);
-        int[] nPtsReq = Strategy.getNPtsReq(dtyp);     // Number of points required for each pattern
-        rTol = Strategy.getRTol(dtyp);  // max. distance from origin in xy
-        zTol = Strategy.getZTol(dtyp);  // max distance from origin in z
-        setupStrategy(MaxLevel,minRadius,Npatterns,triplets,nPtsReq);

 		Detector d = Detector.instance();
 		PropertySet smparm = d.getParameters("HitSmearing");
 	       // ~central tracker resolution in z (in fact length of strip segm).
 		ctrzres = 3.464 * smparm.getDouble("TPCZSmear");

-		vxterr = smparm.getDouble("CCDXYSmear");
-		System.out.println("ctrzres="+ctrzres+" vxdres="+vxterr);
-                

+		vxterr = 3. * smparm.getDouble("CCDXYSmear");
+//		System.out.println("ctrzres="+ctrzres+" vxdres="+vxterr);

 /*		if(dohist)
 		{
 	        String folder = "/Tracking/"+tracker.getType()+"/"+Finder;
 	        HistogramFolder.setDefaultFolder(folder);

-			h010 = new Histogram("Chis/dof of vxd hits");

+			h4 = new Histogram("Dist in Z to corr hit in lr4 vxd");
+			h010 = new Histogram("Chis/dof of all vxd hits");
+			h010_0 = new Histogram("Chis/dof of lr 0 vxd hits");
+			h010_1 = new Histogram("Chis/dof of lr 1 vxd hits");
+			h010_2 = new Histogram("Chis/dof of lr 2 vxd hits");
+			h010_3 = new Histogram("Chis/dof of lr 3 vxd hits");
+			h010_4 = new Histogram("Chis/dof of lr 4 vxd hits");
+			h010_5 = new Histogram("Chis/dof of lr 0 EC vxd hits");
+			h010_6 = new Histogram("Chis/dof of lr 1 EC vxd hits");
+			h010_7 = new Histogram("Chis/dof of lr 2 EC vxd hits");
+			h010_8 = new Histogram("Chis/dof of lr 3 EC vxd hits");

 			h011 = new Histogram("Chis/dof of tracker hits");
 			h012 = new Histogram("Nhits in central tracker for recon. tracks");

+			h012a = new Histogram("Nhits in VXD for recon traks");
+			h012am = new Histogram("Missing lr. nmbrs in VXD for recon traks");
+			h012az = new Histogram("Z-pos. of missing VXD lr 1 hit for recon traks");
+			h012nc = new Histogram("Chi2 of Layer 1 CT hit to VXD trk");
+			h012d = new Histogram("Dist to correct hit in VXD lr 1 when it can't link");
+			h112  = new Histogram2D("N VXD hits vs cos(theta)");
+			h012b = new Histogram("Total nmb of assigned hits for recon traks");

 			h013 = new Histogram("norm.dist with Rcur ~ 100. cm ");
 			h014 = new Histogram("norm.dist with Rcur < 50. cm");
 			h015 = new Histogram("CT residual (microns) in lr1 Rcur > 50 cm");

@@ -274,8 +341,10 @@

 			h022 = new Histogram("assign hits residual in lr2");
 			h023 = new Histogram("assign hits residual in lr3");
 			h024 = new Histogram("assign hits residual in lr4");

+			h036 = new Histogram("Extra VXD layer residual for EC LR 6");
+			h037 = new Histogram("Extra VXD layer residual for EC LR 7");
+			h038 = new Histogram("Extra VXD layer residual for EC LR 8");

 		} */

-

        }
     }
 

@@ -285,6 +354,14 @@

 	}
     public void findTracks(int level, int nPoints,int nLayers)
     {

+        int MaxLevel = Strategy.getLevels(dtyp) - 1;
+        double[] minRadius = Strategy.getMinRadius(dtyp);
+        int Npatterns = Strategy.getNPatterns(dtyp);
+        int[][] triplets = Strategy.getTriplets(dtyp);
+        int[] nPtsReq = Strategy.getNPtsReq(dtyp);     // Number of points required for each pattern
+        rTol = Strategy.getRTol(dtyp);  // max. distance from origin in xy
+        zTol = Strategy.getZTol(dtyp);  // max distance from origin in z
+        setupStrategy(MaxLevel,minRadius,Npatterns,triplets,nPtsReq);

         x=tracker.getX();
         y=tracker.getY();
         z=tracker.getZ();

@@ -296,7 +373,7 @@

         Level = level;
  
 	if (level > MaxLevel) {

-	    System.out.println(Finder+"::findTracks: Max. level is "+MaxLevel);

+	    System.err.println(Finder+"::findTracks: Max. level is "+MaxLevel);

 	    return;
 	}
 	int iMin = 0;

@@ -304,13 +381,13 @@

 	if (level == 0) 
     {
 	 nEvt++;

-     System.out.println("VXDPatFind: MaxPoints= "+MaxPoints);

+//     System.out.println("VXDPatFind: MaxPoints= "+MaxPoints);

 	 mcptot = null;
 	 ntra  = 0;
 	 ntras = 0;
      if(nPoints > MaxPoints)
      {

-      System.out.println("VXDPatFind:: too many points!! "+nPoints+" > "+MaxPoints);

+      System.err.println("VXDPatFind:: too many points!! "+nPoints+" > "+MaxPoints);

       return;
      }
      if(debug) System.out.println("VXDPatFind.findTracks("+level+","+nPoints+","+nLayers+")");

@@ -321,20 +398,30 @@

      // ispin for the given layer/phi bin, and ispon keeps number of entries for the given
      // bin  
      int nlt = nvxlrs+nctlrs+neclrs;

+	 int[] nppll = new int[nlt];
+	 for(int i=0; i<nlt; i++)
+	 {
+	  nppll[i] = 0;	 
+	  for(int j=0; j<100; j++)
+	   nppll[i]+=ilon[i*100+j];
+	 } 

      int nse = nlt*100;
      for(int ind=0; ind<nse; ind++)
       ispon[ind] = 0;

+	 mcphc = -1; 

      for(int lr = 0; lr<nlt; lr++)
      {
 	  int slr = lr*100; 	 

-      int np = ipoi[slr+100]-ipoi[slr];

+//      int np = ipoi[slr+100]-ipoi[slr];
+      int np = nppll[lr];
+	  if((traceAll) && (lr == (nlt-1))) System.out.println("Have "+np+" hits in outer most lr "+lr);

       for(int n=0; n<np; n++)
       { 
        int ip = ipoi[slr]+n;

-	   if(trace && (ip==tracehc))

+	   if(trace && ((ip==tracehc) || traceAll))

 	   {
 	    mcptot = Thits[ip].getMCParticle();

-	    System.out.println("found VXD hit to trace");

+	    System.out.println("found hit to trace");

 		if(mcptot != null) mcphc = mcptot.hashCode();
 	    System.out.println("It's hash code is: "+mcphc);
 	    if(mcptot == null) System.out.println("but it has not defined MC particle");

@@ -355,7 +442,21 @@

        double pphi = Math.atan2(y[ip],x[ip]);
        if(pphi<0.) pphi+=TwoPI;
        int is = (int) ((pphi*100.)/TwoPI);

+	   if((traceAll) && (lr == (nlt-1))) System.out.println("Phi index of last lr hit is "+is);

        if(is < 100) ispon[slr+is]++;

+	   if(trace && (mcphc != -1))
+	   {
+		 MCParticle hmcp = Thits[ip].getMCParticle();
+		 if(hmcp != null)
+		 {
+		  if(hmcp.hashCode() == mcphc)
+		  {
+		   double[] hp=Thits[ip].getPoint();	  
+		   System.out.println("MC with hashcode "+mcphc+" has hit in layer "+ lr+" : "+
+		   df.format(hp[0])+" : "+df.format(hp[1])+" : "+df.format(hp[2])+" phi ind: "+is);
+	      }   
+		 } 
+	   }   

       }
      }
      ispoi[0]=0;

@@ -366,7 +467,8 @@

      for(int lr = 0; lr<nlt; lr++)
      {
 	  int slr = lr*100;	 

-      int np = ipoi[slr+100]-ipoi[slr];

+//      int np = ipoi[slr+100]-ipoi[slr];
+      int np = nppll[lr];

       for(int n=0; n<np; n++)
       { 
        int ip = ipoi[slr]+n;

@@ -401,60 +503,78 @@

 	for (int i = iMin; i < Npatterns; i++) {
             npchkd = 0;
             Paternmb = i; 

-	    int outerLayer  = triplets[i][0];
-	    int middleLayer = triplets[i][1];
-	    int innerLayer  = triplets[i][2];
-		if(trace) System.out.println("At level: "+level+" patern: "+i+" triplet is: "+
-		innerLayer+","+middleLayer+","+outerLayer);
-            extraLayer1 = nvxlrs-1;
-            extraLayer2 = nvxlrs-1;
-            for(int n=0; n<nvxlrs; n++)
-             if((n!=innerLayer)&&(n!=middleLayer)&&(n!=outerLayer)) extraLayer1 = n; 
-            for(int n=0; n<nvxlrs; n++)
-             if((n!=innerLayer)&&(n!=middleLayer)&&(n!=outerLayer)&&(n!=extraLayer1)) extraLayer2 = n; 
-            MaxLayerToCount = extraLayer1; 

+	        outerLayer  = triplets[i][0];
+	        middleLayer = triplets[i][1];
+	        innerLayer  = triplets[i][2];
+//		dotrace = (outerLayer > (nvxblrs-1));
+//		trace = dotrace;
+//		if(dotrace) ntracepr++;
+//		if(i > 1) System.out.println("At level: "+level+" patern: "+i+" triplet is: "+
+//		innerLayer+","+middleLayer+","+outerLayer);
+			if(useVXEC)
+			{
+			 if(outerLayer > nvxblrs)
+			 {
+			 } 
+			}
+            MaxLayerToCount = nvxlrs-1; 

        if(debug) System.out.println("level "+level+" extr1="+extraLayer1+" extr2="+extraLayer2);
 	    // For this given triplet of layers, find all possible tracks
 

-	    selectThreePoints(nPoints, nLayers,
-			      radiusMin, innerLayer, middleLayer, outerLayer);

+	    selectThreePoints(nPoints, nLayers,radiusMin);

             
 	    ntras = ntra;  //prepare ntras for next triplet iteration
 
 	}//end the triplet i loop
 	for(int i=0; i<ntra; i++) fitsuccess[i]=false;     

-        if((level==MaxLevel) && debug)

+        if((level==MaxLevel) && trace)

         { 
          System.out.println("Number of tracks found = "+ntra);
          for(int i=0; i<ntra; i++)
          {

-          System.out.println(" Track "+i+" phi="+df.format(get_phi0(i))+" tan(Lam)="+
-          df.format(get_tandip(i))+" curv="+df.format(get_curv(i))+" n points="+

+          System.out.println(" Track "+i+"trk parms: "+df.format(get_d0(i))+" ; "+
+		  df.format(get_phi0(i))+" ; "+df.format(get_curv(i))+" ; "+
+		  df.format(get_z0(i))+" ; "+df.format(get_tandip(i))+" n points="+

           getNPoints(i));

+		  System.out.print("pnts: ");
+		  for(int j=0; j<getNPoints(i); j++)
+		   System.out.print(itra[i][j+1]+" ");
+		  System.out.println(""); 

          }
         }
 	return;
     }
 

-    protected void selectThreePoints(int nPoints,int nLayers,double radiusMin, int innerLayer, 
-                                   int middleLayer, int outerLayer)

+    protected void selectThreePoints(int nPoints,int nLayers,double radiusMin)

     {
         boolean first = false;

+		boolean innerEC = (innerLayer > (nvxblrs - 1));
+		boolean outerEC = (outerLayer > (nvxblrs - 1));
+		boolean middleEC = (middleLayer > (nvxblrs - 1));

 	/** Selects three suitable points for a given triplet of layers */
         int ntracan = 0;

-        int[] vxdhin = new int[nvxlrs];
-		int[] ntamb = {3,2,1,1,1};

+		int[] ntamb = {10,7,5,3,3,5,5,5,5};
+		double[] vxchisq= new double[9];

 
         if((innerLayer > nLayers)||(outerLayer > nLayers) || (middleLayer > nLayers) ||
           (innerLayer < -1) || (outerLayer < 1) || (middleLayer < 0))
         {

-         System.out.println("selectThreePoints:Error in Layer number");

+         System.err.println("selectThreePoints:Error in Layer number.We have only "+nLayers);

          return;
         }

+        if((innerLayer >=nvxlrs)||(outerLayer >=nvxlrs) || (middleLayer >= nvxlrs))
+         return;
+		if((!useVXEC)&&(innerEC || middleEC || outerEC)) return; 

 	int howMany=0;
 	int istrt = innerLayer*100;
 	int iend = istrt+100;
 	int nhinn = 1;

+	for(int i=0; i<5; i++)
+	 ntamb[i]= (int) (10. * (zTol/vdr[i]));
+	ntamb[0]+=2;
+	ntamb[1]+=1;
+	if(debug) System.out.println("Tnlam bins: "+ntamb[0]+" "+ntamb[1]+" "+ntamb[2]+" "+
+	 ntamb[3]+" "+ntamb[4]);

 	if (innerLayer != -1) 
 	{
 	  nhinn = ipoi[iend]-ipoi[istrt];

@@ -466,6 +586,8 @@

 	iend=istrt+100;
 	int nhout = ipoi[iend]-ipoi[istrt];
 	howMany = nhinn*nhmidd*nhout;

+	if(trace) System.out.println("Layers: "+innerLayer+" "+middleLayer+" "+outerLayer); 
+	if(trace) System.out.println("Have "+nhinn+" inner hits, "+nhmidd+" middle "+nhout+" outer");

 /**/
 	// *   Loop over sets of triplets of layers,
 	// * only if there are points in ALL layers. 

@@ -495,6 +617,11 @@

 	  x[initialInner] = y[initialInner] = z[initialInner] = 0.;
 	  KFlag[initialInner] = false;
 	 }

+	 if(trace)
+	 {
+	  for(int i=initialInner; i<finalInner; i++)
+       if(Thits[i].getMCParticle() == mcptot) System.out.println("inner layer "+innerLayer+" HAS traced hit");	  
+	 } 

      for(int ntbo=0; ntbo<100; ntbo++)
 	 {
 	  initialOuter = ipoi[outerLayer*100+ntbo];

@@ -520,38 +647,66 @@

 	  double leverarmrp = 2.4;
 	  double leverarmf = 2.4;
 	  for (int IOuter = initialOuter; IOuter <= finalOuter; IOuter++) 

-      {              

+      {
+       if(debug) System.out.println("Outer hit index: "+IOuter);
+       if(!trace) mcptot = Thits[IOuter].getMCParticle();              

        //loop over points in inner layer
 	   for (int IInner = initialInner; IInner <= finalInner; IInner++) 
        {

-		if(Math.abs(MinDPhi(mainco[IOuter],mainco[IInner]))<eighthPI)

+		if(debug) System.out.println("Inner hit index: "+IInner);   
+//		if((Math.abs(MinDPhi(mainco[IOuter],mainco[IInner]))<eighthPI) || outerEC)
+        if(true)

         { 
          //loop over points in middle layer
          for (int IMiddle = initialMiddle; IMiddle <= finalMiddle; IMiddle++) 
          {

-	      if(Math.abs(MinDPhi(mainco[IOuter],mainco[IMiddle])) < eighthPI)

+		  sameMCPhts = false;
+		  fromIP = false;
+		  if(debug) System.out.println("Middle hit index: "+IMiddle);	 
+//	      if((Math.abs(MinDPhi(mainco[IOuter],mainco[IMiddle])) < eighthPI) || outerEC)
+          if(true)

           {
            // First (soft) z-check - middle point z is between inner z and outer z
 		   dotrace = false;

-           if((trace) && (mcptot != null))

+           if(mcptot != null)

            {

+			MCParticle mcpinr =  Thits[IInner].getMCParticle();
+			MCParticle mcpmdl = Thits[IMiddle].getMCParticle();
+			MCParticle mcpout = Thits[IOuter].getMCParticle();
+			int mcphcout = 0;
+			int mcphcinr = 0;
+			int mcphcmdl = 0;
+			if(mcpout != null) mcphcout = mcpout.hashCode();
+			if(mcpinr != null) mcphcinr = mcpinr.hashCode();
+			if(mcpmdl != null) mcphcmdl = mcpmdl.hashCode();
+//			System.out.println("hash codes of hits MCP: "+mcphcinr+" "+mcphcmdl+" "+mcphc);

 		    if((Thits[IOuter].getMCParticle()==mcptot) && 
 			 (Thits[IInner].getMCParticle() == mcptot) &&

-			  (Thits[IMiddle].getMCParticle() == mcptot)) dotrace=true;
-		   } 

+			  (Thits[IMiddle].getMCParticle() == mcptot)) sameMCPhts=true;
+			if(sameMCPhts)
+			{
+			 double[] coo = mcptot.getOrigin();
+			 if(Math.sqrt(coo[0]*coo[0]+coo[1]*coo[1]+coo[2]*coo[2]) < 0.01) fromIP = true;
+			}
+			int ncoin = 0;
+			if(mcphc == mcphcout) ncoin++;
+			if(mcphc == mcphcinr) ncoin++;
+			if(mcphc == mcphcmdl) ncoin++;
+			if((ncoin > 2) && (trace))
+			 System.out.println("("+innerLayer+" "+middleLayer+" "+outerLayer+") inner "+mcphcinr+" middle "+mcphcmdl+" outer "+mcphcout);
+			if(trace && sameMCPhts) dotrace=true;
+		   }
+		   if((traceAll) && sameMCPhts) dotrace=true;

 		   if(dotrace) System.out.println("tracing reconstruction: iout= "+IOuter+" iinn= "+
 			IInner+" imid= "+IMiddle);
            boolean zCheck = false;
            if(z[IInner] <= z[IOuter])

-            if((z[IMiddle] >= z[IInner])&&(z[IMiddle]<=z[IOuter])) zCheck = true;

+            if((z[IMiddle] >= z[IInner]-0.1)&&(z[IMiddle]<=z[IOuter]+0.1)) zCheck = true;

            if(z[IInner] >= z[IOuter])

-            if((z[IMiddle] >= z[IOuter])&&(z[IMiddle] <= z[IInner])) zCheck = true; 

+            if((z[IMiddle] >= z[IOuter]-0.1)&&(z[IMiddle] <= z[IInner]+0.1)) zCheck = true;
+		   if(dotrace) System.out.println("Z CHECK : "+zCheck);	

            if(zCheck)
            {

-		    leverarmrp = Math.sqrt((x[IOuter]-x[IInner])*(x[IOuter]-x[IInner]) +
-			 (y[IOuter]-y[IInner])*(y[IOuter]-y[IInner]));
-			leverarmf = Math.sqrt(leverarmrp*leverarmrp+
-			 (z[IOuter]-z[IInner])*(z[IOuter]-z[IInner]));

 			if(dotrace) System.out.println("Pass first zCheck");	  
             p1[0]=x[IInner];
             p1[1]=y[IInner];

@@ -560,11 +715,17 @@

             p3[0]=x[IOuter];
             p3[1]=y[IOuter];
             find3ptrack(p1,p2,p3,z[IInner],z[IOuter],true);

+			if(dotrace) System.out.println("R curvature of trk.cand: "+df.format(Rcu)+" while minimum is: "+
+			 radiusMin);

 		    boolean curvCheck = (Rcu > radiusMin);
 		    if(curvCheck) 
             {

+		     leverarmrp = 0.5 * Math.sqrt((x[IOuter]-x[IInner])*(x[IOuter]-x[IInner]) +
+			 (y[IOuter]-y[IInner])*(y[IOuter]-y[IInner]));
+			 leverarmf = Math.sqrt(leverarmrp*leverarmrp+
+			  0.25 * (z[IOuter]-z[IInner])*(z[IOuter]-z[IInner]));

 		     if(dotrace) System.out.println("pass curv. check");

-			 double smlever = 0.6;

+			 double smlever = 1.2;

 			 double angtersq = (vxterr*vxterr)/(leverarmrp*leverarmrp);
 			 double mstermsq = vxdmsco*vxdmsco*trkprms.omega*trkprms.omega;
 			 double angerr = Math.sqrt(angtersq + mstermsq);

@@ -574,23 +735,33 @@

 			 vxdmsco*vxdmsco*trkprms.omega*trkprms.omega);
 			 vxrpersq = angerr*angerr*smlever*smlever+vxterr*vxterr;
 			 vxdphersq = angerr*angerr;

-			 double zersc = (1.+Math.abs(trkprms.s))*smlever*angerz;

+//			 double zersc = (1.+Math.abs(trkprms.s))*smlever*angerz;
+			 double zersc = (1.+trkprms.s*trkprms.s)*(1.+trkprms.s*trkprms.s)*smlever*angerz;

 			 vxtolsq = zersc*zersc+vxterr*vxterr;
 			 zerv = Math.sqrt(vxtolsq);
 			 zerv = 5.*zerv;
              // now we can do more precise z-check for middle point
              double phii = Math.atan2(p1[1]-yc[ntra],p1[0]-xc[ntra]);
              if(phii < 0.) phii += TwoPI;

+			 if(phii > TwoPI) phii -= TwoPI;

              double phim = Math.atan2(p2[1]-yc[ntra],p2[0]-xc[ntra]);
              if(phim < 0.) phim += TwoPI;

+			 if(phim > TwoPI) phim -= TwoPI;

              double phio = Math.atan2(p3[1]-yc[ntra],p3[0]-xc[ntra]); 
              if(phio < 0.) phio += TwoPI;

+			 if(phio > TwoPI) phio-= TwoPI;

              double dph1 = Math.abs(phim-phii);
              if(dph1 > Math.PI) dph1 = Math.abs(dph1-TwoPI);
              double dph2 = Math.abs(phio-phim);
              if(dph2 > Math.PI) dph2 = Math.abs(dph2-TwoPI);
              double zexp = z[IInner] + (z[IOuter]-z[IInner]) * dph1/(dph1+dph2);
              double zmis = Math.abs(z[IMiddle]-zexp);

+			 if(zmis > zerv)
+			 {
+			  dph2 = TwoPI - dph2;	 
+              zexp = z[IInner] + (z[IOuter]-z[IInner]) * dph1/(dph1+dph2);
+              zmis = Math.abs(z[IMiddle]-zexp);
+			 } 

              if(zmis < zerv) //middle point z  close to expected z
              {
 			  if(dotrace) System.out.println("pass tight Z check");	

@@ -606,84 +777,214 @@

 				lstvhphi = mainco[IOuter];
 			    if(dotrace) System.out.println("pass d0 check");	  
                 swmr.swim(trkprms);

-                // we will look now for hits on this track in 2 more VXD layers,

+                // we will look now for hits on this track in  VXD layers

                 // not included in the triplet

-                int nexthts = 0;

+                nexthts = 0;

                 boolean havhit = false;

-                int[] exlrs = new int[2];
-                exlrs[0]=extraLayer1;
-                exlrs[1]=extraLayer2;

 			    ndegfr = 0;
 			    trvxchsq = 0.;

-                for(int extr=0; extr<2; extr++)

+//				System.out.println("Looking for extra hits in VXD");
+                double zpozm=0.;
+				double dtoch=-1.;
+				double xpozm=0.;
+				double ypozm = 0.;
+				double anybd = 10.;
+				int crL1hi = -1;
+				boolean trhtl1 = false;
+                for(int lrn=0; lrn<nvxlrs; lrn++)

                 {

-                 double bestdis = 25.;
-                 havhit = false;
-                 int lrn = exlrs[extr];
-                 vxdhin[lrn]=-1;
-                 double rts = vspecs.getLayerRadius(lrn);
-                 double zts = vspecs.getLayerZ(lrn);
-                 swmr.swim(rts,zts);
-                 int plane = swmr.getPlane();
-                 if(plane == swmr.PLANE_R)
-                 {
-                  double[] xpn = swmr.getIntersect();
-                  double phpn = Math.atan2(xpn[1],xpn[0]);
-                  if(phpn < 0.) phpn+=TwoPI;
-                  int iph = (int) ((phpn*100.)/TwoPI);
-                  int[] isr = new int[3];
-                  isr[0] = iph-1;
-                  if(isr[0] < 0) isr[0] = 99;
-                  isr[2] = iph+1;
-                  if(isr[2] > 99) isr[2] = 0;
-                  isr[1] = iph;
-                  for(int ns=0; ns<3; ns++)
-                  {
-                   int is = lrn*100+isr[ns];
-                   int np = ispon[is];
-                   for(int pn=0; pn<np; pn++)

+				 vxchisq[lrn]=-1.;	
+				 if((lrn!=outerLayer)&&(lrn!=middleLayer)&&(lrn!=innerLayer))
+				 {
+				  int lvri = Math.abs(innerLayer-lrn);
+				  if(lvri > 5) lvri-=5;
+				  int lvrm = Math.abs(middleLayer-lrn);
+				  if(lvrm > 5) lvrm-=5;
+				  int lvro = Math.abs(outerLayer-lrn);
+				  if(lvro > 5) lvro-=5;
+				  int lvrco = lvri;
+				  if(lvrm < lvrco) lvrco = lvrm;
+				  if(lvro < lvrco) lvrco = lvro;
+				  double lvrfa = lvrco;
+				  if((lrn > innerLayer) && (lrn < outerLayer)) lvrfa = 0.5 * lvrco;
+				  if(lvrfa < 0.5) System.out.println ("ERROR in determining LEVER factor");
+                  double bestdis = 10.;
+				  if(lrn == 0) anybd = 10.;
+                  havhit = false;
+                  vxdhin[lrn]=-1;
+				  if(lrn < nvxblrs)
+				  {
+                   double rts = vspecs.getLayerRadius(lrn);
+                   double zts = vspecs.getLayerZ(lrn);
+                   swmr.swim(rts,zts);
+                   int plane = swmr.getPlane();
+                   if(plane == swmr.PLANE_R)

                    {

-                    int ip = ispoi[is]+pn;
-                    int ii = ispin[ip];
-                    double xde = xpn[0]-x[ii];
-                    double yde = xpn[1]-y[ii];
-				    double xych2 = (xde*xde+yde*yde)/vxrpersq;
-                    double zde = xpn[2]-z[ii];
-				    double zch2 = zde*zde/vxtolsq;
-                    double dds = xych2+zch2;
-                    if(dds<bestdis)
-                    {  
-                     havhit=true;
-                     bestdis=dds;
-                     vxdhin[lrn]=ii;
-                    }   

+                    double[] xpn = swmr.getIntersect();
+                    double phpn = Math.atan2(xpn[1],xpn[0]);
+					if(lrn == 0)
+					{
+					 trhtl1 = true;	
+					 xpozm = xpn[0];
+					 ypozm = xpn[1];
+					 zpozm = xpn[2];
+				    } 
+                    if(phpn < 0.) phpn+=TwoPI;
+                    int iph = (int) ((phpn*100.)/TwoPI);
+                    int[] isr = new int[5];
+                    isr[0] = iph-2;
+                    if(isr[0] < -1) isr[0] = 98;
+					if(isr[0] < 0) isr[0]= 99;
+					isr[1]=iph-1;
+					if(isr[1] < 0) isr[1] = 99;
+                    isr[3] = iph+1;
+                    if(isr[3] > 99) isr[3] = 0;
+					isr[4] = iph+2;
+					if(isr[4] > 100) isr[4] = 1;
+					if(isr[4] > 99) isr[4] = 0; 
+                    isr[2] = iph;
+                    for(int ns=0; ns<5; ns++)
+                    {
+                     int is = lrn*100+isr[ns];
+                     int np = ispon[is];
+                     for(int pn=0; pn<np; pn++)
+                     {
+                      int ip = ispoi[is]+pn;
+                      int ii = ispin[ip];
+ 					  MCParticle hmc = Thits[ii].getMCParticle();
+                      double xde = xpn[0]-x[ii];
+                      double yde = xpn[1]-y[ii];
+				      double xych2 = (xde*xde+yde*yde)/(vxrpersq*lvrfa*lvrfa);
+                      double zde = xpn[2]-z[ii];
+				      double zch2 = zde*zde/(vxtolsq*lvrfa*lvrfa);
+                      double dds = xych2+zch2;
+/*					  if((sameMCPhts) && (hmc==mcptot) && (lrn == 4)&& (Math.abs(trkprms.s) >0.9) &&
+					   (Math.abs(trkprms.s) < 1.1))
+					   h4.fill(zde); */
+					  if(sameMCPhts && (lrn == 0))
+					  {
+					   double ahdist = Math.sqrt(xde*xde+yde*yde+zde*zde);
+					   if(ahdist < anybd) anybd = ahdist;
+					  }  
+					  if(sameMCPhts && (lrn==0) && (hmc == mcptot))
+					  {
+					    dtoch = Math.sqrt(dds);
+						crL1hi = ii;
+					  }
+					  if(dotrace && (hmc==mcptot))
+					  {
+					   System.out.println("VXD lr: "+lrn+" nrm dst.sq to corr hit: "+df.format(dds)+
+					   " abs dis in X: "+df.format(xde)+" in Y: "+df.format(yde)+" in Z: "+df.format(zde));
+					   System.out.println(" exp. X: "+df.format(xpn[0])+" Y: "+df.format(xpn[1])+
+					   " Z: "+df.format(xpn[2]));
+				      }   
+                      if(dds<bestdis)
+                      {  
+                       havhit=true;
+                       bestdis=dds;
+                       vxdhin[lrn]=ii;
+                      }
+				     }  
+                    }

                    }
                   }

-                 }
-                 if(havhit)
-				 { 
-				  nexthts++;
-				  trvxchsq+=bestdis;
-				  }
-                 } 
-                 if(nexthts >= nextvh)

+				  else
+				  {
+				   if(useVXEC)
+				   {
+				    int elrn = lrn-nvxblrs;	  
+                    double rts = vespecs.getLayerOuterRadius(elrn);
+				    double rmax = specs.getOuterRadius();
+                    double zts = vespecs.getLayerZ(elrn);
+                    swmr.swim(rmax,zts);
+                    int plane = swmr.getPlane();
+				    double evxrpersq = 0.025 * trkprms.omega * trkprms.omega*(1.+trkprms.s*trkprms.s);
+                    if(((plane==swmr.PLANE_PLUSZ) || (plane==swmr.PLANE_MINUSZ))&&useVXEC)
+				    {
+                     if(dotrace) System.out.println("Helix hit VXD disk! ");
+					 if(dotrace) System.out.println("VX error sq: "+evxrpersq);
+					 double[] xpn = swmr.getIntersect();
+					 double her = Math.sqrt(xpn[0]*xpn[0]+xpn[1]*xpn[1]);
+                     double phpn = Math.atan2(xpn[1],xpn[0]);
+					 if(her < rts)
+					 {
+                      if(phpn < 0.) phpn+=TwoPI;
+                      int iph = (int) ((phpn*100.)/TwoPI);
+                      int[] isr = new int[3];
+                      isr[0] = iph-1;
+                      if(isr[0] < 0) isr[0] = 99;
+                      isr[2] = iph+1;
+                      if(isr[2] > 99) isr[2] = 0;
+                      isr[1] = iph;
+                      bestdis = 10.;
+					  double bestcorhd =10000.;
+                      for(int ns=0; ns<3; ns++)
+                      {
+                       int is = lrn*100+isr[ns];
+                       int np = ispon[is];
+                       for(int pn=0; pn<np; pn++)
+                       {
+                        int ip = ispoi[is]+pn;
+                        int ii = ispin[ip];
+                        double xde = xpn[0]-x[ii];
+                        double yde = xpn[1]-y[ii];
+				        double dds = (xde*xde+yde*yde)/(evxrpersq*lvrfa*lvrfa);
+					    if(sameMCPhts)
+					    {
+						 if(Thits[ii].getMCParticle() == mcptot)
+						 {
+						  double dist = Math.sqrt(dds);
+						  if(dist < Math.abs(bestcorhd)) bestcorhd = dist;
+						  if(xde < 0.) bestcorhd = -bestcorhd;
+						 }
+					    }   
+                        if(dds<bestdis)
+                        {  
+                         havhit=true;
+                         bestdis=dds;
+                         vxdhin[lrn]=ii;
+                        }
+				       }  
+                      }
+/*					  if(Math.abs(bestcorhd) < 10000.)
+					  {
+					   if(elrn == 1) h036.fill(bestcorhd);	 
+					   if(elrn == 2) h037.fill(bestcorhd);	 
+					   if(elrn == 3) h038.fill(bestcorhd);	 
+					  } */
+				     }
+				    }
+				    if((dotrace) && (havhit)) System.out.println("VXDEC: Norm dsq="+bestdis);
+				   }
+			      }   
+                  if(havhit)
+				  {
+				   nexthts++;
+				   trvxchsq+=bestdis;
+				   vxchisq[lrn]=bestdis/2.;
+                  }
+			     } 
+			    }
+                if(nexthts >= nextvh)
+                {
+				 ndegfr=2*nexthts;
+//				 ndegfrvtx = 2*ndegfr;
+                 ndegfrvtx = ndegfr;
+				 if(dotrace) System.out.println("pass extra hit test");   
+                 vxdhin[innerLayer]=IInner;
+                 vxdhin[middleLayer]=IMiddle;
+                 vxdhin[outerLayer]=IOuter;
+                 itra[ntra][0]=0;
+                 for(int nvl=0; nvl<nvxlrs; nvl++)

                  {

-				  ndegfr=2*nexthts;
-				  ndegfrvtx = 2*ndegfr;
-				  if(dotrace) System.out.println("pass extra hit test");   
-                  vxdhin[innerLayer]=IInner;
-                  vxdhin[middleLayer]=IMiddle;
-                  vxdhin[outerLayer]=IOuter;
-                  itra[ntra][0]=0;
-                  for(int nvl=0; nvl<nvxlrs; nvl++)

+                  if(vxdhin[nvl] !=-1)

                   {

-                   if(vxdhin[nvl] !=-1)
-                   {
-                    itra[ntra][0]++;
-                    int inh=itra[ntra][0];
-                    itra[ntra][inh]=vxdhin[nvl];
-                   }
-                  }  

+                   itra[ntra][0]++;
+                   int inh=itra[ntra][0];
+                   itra[ntra][inh]=vxdhin[nvl];
+				   lstvxhin = vxdhin[nvl];
+                  }
+                 }  

 	              // We have a track - increment track counter.
                   ntracan++;
 	              ntra++;

@@ -705,18 +1006,80 @@

 				    if(dotrace) System.out.println("But duplicate!");	 
                     ntra--;
                    } // duplicate check;

+				   if(dotrace)
+				   {
+				    System.out.println("have tracks: "+ntra);
+					System.out.println("track list as it looks now: ");
+					for(int n=0; n<ntra; n++)
+					{
+					 int nntp = itra[n][0];	
+					 System.out.print("Tr: "+n+" npnts "+nntp+" : ");
+					 for(int k=0; k<nntp; k++)
+					  System.out.print(itra[n][k+1]+" ");
+					 System.out.println(""); 
+					}
+				   }

 				   if(ntrat==ntra)
 				   {
 					double efomega = Math.abs(trkprms.omega);
 					if(efomega < 0.01) efomega = 0.01;

-                                        
-/*					if(dohist)

+					if((dohist) && sameMCPhts && fromIP)

 					{
 //					 h101.fill(1./efomega,trvxchsq/ndegfrvtx);

-					 h010.fill(trvxchsq/ndegfrvtx);

+/*					 h010.fill(trvxchsq/ndegfrvtx);
+					 if(vxchisq[0]>0.) h010_0.fill(vxchisq[0]);
+					 if(vxchisq[1]>0.) h010_1.fill(vxchisq[1]);
+					 if(vxchisq[2]>0.) h010_2.fill(vxchisq[2]);
+					 if(vxchisq[3]>0.) h010_3.fill(vxchisq[3]);
+					 if(vxchisq[4]>0.) h010_4.fill(vxchisq[4]);
+					 if(vxchisq[5]>0.) h010_5.fill(vxchisq[5]);
+					 if(vxchisq[6]>0.) h010_6.fill(vxchisq[6]);
+					 if(vxchisq[7]>0.) h010_7.fill(vxchisq[7]);
+					 if(vxchisq[8]>0.) h010_8.fill(vxchisq[8]);

 	                 if(nastrhits > 0) h011.fill(asshchsq/(double) (nastrhits));
 					// if((Math.abs(trkprms.omega) < 0.008)&&(Math.abs(trkprms.s)<1.2))

-					 h012.fill(nastrhits+1);

+					 h012.fill(nastrhits);
+					 h012a.fill(nexthts+3);
+					 for(int n=0; n<nvxlrs; n++)
+					  if(vxdhin[n] == -1) h012am.fill(n+1); */
+					 if(vxdhin[0] == -1)
+					 {
+					  if(sameMCPhts && trhtl1)
+					  {
+				//	   h012az.fill(zpozm);
+					   double[] opn = Thits[IOuter].getPoint();
+					   double[] mpn = Thits[IMiddle].getPoint();
+					   double[] ipn = Thits[IInner].getPoint();
+					   double[] mcmo = mcptot.getMomentum();
+					   double mcPt = Math.sqrt(mcmo[0]*mcmo[0]+mcmo[1]*mcmo[1]);
+					   double mcTl = mcmo[2]/mcPt;
+/*					   System.out.println("Reconstracting real track for MC: "
+					   +mcptot.getType().getName()+" origin: X="+mcptot.getOrigin()[0]+" Y= "+
+					   mcptot.getOrigin()[1]+" Z= "+mcptot.getOrigin()[2]+

[truncated at 1000 lines; 1135 more skipped]

lcd/hep/lcd/recon/tracking/ccd

diff -u -r1.2 -r1.3
--- VXDReco.java	15 Mar 2005 00:46:56 -0000	1.2
+++ VXDReco.java	10 Aug 2005 20:28:16 -0000	1.3
@@ -17,6 +17,7 @@

 import hep.lcd.recon.tracking.combined.*;
 import hep.lcd.mc.smear.*;
 import hep.lcd.recon.tracking.ccd.CCDSpecifications;

+import hep.lcd.recon.tracking.ccd.CCDECSpecifications;

 import hep.lcd.recon.tracking.tpc.TPCSpecifications;
 import hep.lcd.recon.tracking.tpc.EndCpSpecifications;
 

@@ -37,20 +38,26 @@

 {
     Specifications vspecs;
     Specifications ecspecs;

+	CCDECSpecifications vespecs;

 	VertexDetector VD;

-	public double MinVxdHitSep = 0.001;

+	public double MinVxdHitSep = 0.0;

 	public double VxdHitRerr = 0.2;
 	public double TrkrHitRerr = 0.5;
 	public double EcHitZerr = 0.5;

-	boolean includeEC = false;
-	boolean includeVXDEC = false;
-    int VXDLayers = 5;
-    int CTRLayers = 5;
-    int ECTLayers = 10;

+	public boolean includeEC = true;
+	public boolean includeVXDEC = true;
+	boolean uincVXDEC = true;
+	boolean trace = false;
+	double[] trhco = new double[3];
+    int VXDLayers = 0;
+	int VXDECLayers = 0;
+    int CTRLayers = 0;
+    int ECTLayers = 0;

 	int MaxTracks = 200;
 	int nPoints = 0;

-	int[] nvxpnt; 
-	Histogram h011 = null;

+	int[] nvxpnt;
+	MCParticle mcptt = null;
+/*	Histogram h011 = null;

 	Histogram h012 = null;
 	Histogram h013 = null;
 	Histogram h014 = null;

@@ -60,7 +67,12 @@

 	Histogram h023 = null;
 	Histogram h024 = null;
 	Histogram h025 = null;

+	Histogram2D h100 = null;
+	Histogram2D h101 = null;
+	Histogram2D h200 = null;
+	Histogram2D h201 = null; */

 	boolean dohist = false;

+    DecimalFormat df = new DecimalFormat();

     
     VXDPatFind finder;
 

@@ -71,7 +83,7 @@

     }
     public VXDReco(boolean hist)
     {

-   	this(false,null);

+   	this(hist,null);

     }
     /**  Construct VXD reconstructor (MaxPoints=60000)
      * 

@@ -80,8 +92,10 @@

     {
      MaxPoints = 60000;
      finder = new VXDPatFind(hist,strat);

+	 finder.setTraceAll(false);

      MaxTracks = finder.MaxTracks;
 	 dohist = hist;

+     df.setMaximumFractionDigits(5);

      super.finder=finder;
      x = new double[MaxPoints+1];
      y = new double[MaxPoints+1];

@@ -91,11 +105,12 @@

 
         vspecs = new CCDSpecifications();
         specs = new TPCSpecifications(); // it's just central tracker, not necessary TPC 

-        ecspecs = new EndCpSpecifications(); 

+        ecspecs = new EndCpSpecifications();
+		vespecs = new CCDECSpecifications();

         setupParameters();
         nvxpnt = new int[MaxTracks];
 	    VD = new VXDRecoVD(this);

- 		if(dohist)

+ /*		if(dohist)

 		{
 	        String folder = "/Tracking/VXDReco";
 	        HistogramFolder.setDefaultFolder(folder);

@@ -109,12 +124,27 @@

 			h023 = new Histogram("Nhits in 3-rd CT lr");
 			h024 = new Histogram("Nhits in 4-th CT lr");
 			h025 = new Histogram("Nhits in 5-th CT lr");

-	  }

+			h100 = new Histogram2D("XY of traced hits in VXD");
+			h101 = new Histogram2D("RZ of traced hits in VXD");
+			h200 = new Histogram2D("XY of traced hits in Tracker");
+			h201 = new Histogram2D("RZ of traced hits in Tracker");
+	  } */

    }
 
     /** Get type of track hits. */
     public String getType() { return "COMB"; }
 	

+	public void setTrace(boolean set, double[] pnt)
+	{
+		trace = set;
+		if(trace)
+		{
+			trhco[0]=pnt[0];
+			trhco[1]=pnt[1];
+			trhco[2]=pnt[2];
+		}
+	}
+	

 	public VertexDetector getVertexDetector()
 	{
 		return VD;

@@ -128,6 +158,8 @@

     public Specifications getVxdSpecifications() { return vspecs; }
     
     public Specifications getECSpecifications() { return ecspecs; }

+	
+	public Specifications getVxdECSpecifications() { return vespecs; }

   
     public int getMaxPoints() { return MaxPoints; } 
 

@@ -154,18 +186,22 @@

 
 	public void setIncludeVXDEC(boolean include)
 	{

-		includeVXDEC = include;

+		uincVXDEC = include;

 	}
 	public void setupParameters()
 	{

-		finder.setupFinder(this);

                 VXDLayers = vspecs.getNLayers();
                 CTRLayers = specs.getNLayers();
                 ECTLayers = ecspecs.getNLayers();

-	        NLayers  = VXDLayers+CTRLayers+ECTLayers;

+				VXDECLayers = vespecs.getNLayers();
+		includeVXDEC = uincVXDEC && vespecs.haveEndCap();	
+	    if(!includeVXDEC) VXDECLayers = 0;
+		NLayers  = VXDLayers+VXDECLayers+CTRLayers+ECTLayers;
+		finder.setupFinder(this);

 		MaxLevel = finder.getMaxLevel();

-		System.out.println("Total number of layers - VXD: "+VXDLayers+" CT: "+CTRLayers+
-		" ECT: "+ECTLayers);

+		System.out.println("Tot. nmbr of lrs - VXD: "+VXDLayers+" VXEC:"+
+		VXDECLayers+" CT: "+CTRLayers+" ECT: "+ECTLayers);
+		System.out.println("finder MaxLevel is: "+MaxLevel);

 	}
 
 	public int getNPoints(int itrk)

@@ -276,6 +312,7 @@

     {
      specs.update(event);
      vspecs.update(event);

+	 vespecs.update(event);

      ecspecs.update(event);
      if(fitter != null) fitter.setupEvent(event);
      if(specs.areNew()) setupParameters();

@@ -289,8 +326,7 @@

 		}
        VD.startEvent(event);         
 		TrackerHits trkhits = event.getTrackerHits();

-                TrackerHits vxdhits = event.getVXDHits();
-  

+        TrackerHits vxdhits = event.getVXDHits();

 		if (Hist) {
 			String folder = "/Tracking/"+getType();
 			HistogramFolder.setDefaultFolder(folder);

@@ -305,47 +341,107 @@

 		prepnt[1]=0.;
 		prepnt[2]=0.;
 		for (Enumeration hits = vxdhits.getHits();(hits.hasMoreElements() && ncl<MaxPoints);) 

-           {
-			TrackerHit hit = (TrackerHit) hits.nextElement();
-			int layer = hit.getLayer();
-            if(layer > 127) layer-=128;
-			int slayer = layer*100;
-            if(layer < VXDLayers)
-            {
-			 if (Hhits) histogram("Hits/layer").fill(layer);
-            }
-            double specr=vspecs.getLayerRadius(layer);
-			double point[] = hit.getPoint();
-            double acrad=Math.sqrt(point[0]*point[0]+point[1]*point[1]);
-			double dpp = Math.sqrt((point[0]-prepnt[0])*(point[0]-prepnt[0]) +

+        {
+		 TrackerHit hit = (TrackerHit) hits.nextElement();
+		 int layer = hit.getLayer();
+         if(layer > 127) layer-=128;
+		 int slayer = layer*100;
+         if(layer < VXDLayers)
+         {
+		   if (Hhits) histogram("Hits/layer").fill(layer);
+         }
+		 if(!hit.isEndcap())
+		 {
+          double specr=vspecs.getLayerRadius(layer);
+		  double point[] = hit.getPoint();
+		  if(trace)
+		  {
+			if(Math.abs(point[0]-trhco[0]) < 0.0001)
+			{
+			 if(Math.abs(point[1]-trhco[1]) < 0.0001)
+			 {
+			  if(Math.abs(point[2]-trhco[2]) < 0.0001)
+			  {
+				mcptt = hit.getMCParticle();  
+			  }  
+			 } 
+			}
+		  }  
+          double acrad=Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+		  double dpp = Math.sqrt((point[0]-prepnt[0])*(point[0]-prepnt[0]) +

 			(point[1]-prepnt[1])*(point[1]-prepnt[1])+(point[2]-prepnt[2])*(point[2]-prepnt[2]));

-			prepnt[0]=point[0];
-			prepnt[1]=point[1];
-			prepnt[2]=point[2];
-             if(Math.abs(specr-acrad) > VxdHitRerr)

+		  prepnt[0]=point[0];
+		  prepnt[1]=point[1];
+		  prepnt[2]=point[2];
+          if(Math.abs(specr-acrad) > VxdHitRerr)
+		  {
+			 System.out.println("VXD hit position R is: "+acrad+" while spec R: "+specr+
+			 " for layer "+layer);
+		  }
+		  if(dpp < MinVxdHitSep) System.out.println("Too close hits in VXD");
+//          if((Math.abs(specr-acrad) < VxdHitRerr)&&(dpp>MinVxdHitSep))
+          if(Math.abs(specr-acrad) < VxdHitRerr)
+          {
+		   double tle = point[2]/specr;
+		   tle+=5.0;
+		   if(tle<0.) tle=0.;
+		   int tind = (int) (tle*10.);
+		   if(tind>99) tind=99;
+		   tind+=slayer;
+		   nlayer[tind]++;
+		   ncl++;
+          }
+		  else
+		  {
+			System.out.println("VXD hit actual rad: "+acrad+" while spec is: "+specr);  
+		  }  
+		 }
+		 if((hit.isEndcap()) && includeVXDEC)
+		 {
+		  int tlayer = layer;	 
+		  layer+= VXDLayers;
+		  slayer = layer*100;
+		  double point[] = hit.getPoint();
+		  if(trace)
+		  {
+			if(Math.abs(point[0]-trhco[0]) < 0.0001)
+			{
+			 if(Math.abs(point[1]-trhco[1]) < 0.0001)

 			 {

-				 System.out.println("VXD hit position R is: "+acrad+" while spec R: "+specr+
-				 " for layer "+layer);

+			  if(Math.abs(point[2]-trhco[2]) < 0.0001)
+			  {
+				mcptt = hit.getMCParticle();  
+			  }  

 			 } 

-            if((Math.abs(specr-acrad) < VxdHitRerr)&&(dpp>MinVxdHitSep))
-            {
-			 double tle = point[2]/specr;
-			 tle+=5.0;
-			 if(tle<0.) tle=0.;
-			 int tind = (int) (tle*10.);
-			 if(tind>99) tind=99;
-			 tind+=slayer;
-			 nlayer[tind]++;
-			 ncl++;
-            }
-		}

+			}
+		  }  
+          double acrad=Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+          double specz=vespecs.getLayerZ(tlayer);
+          double abz=Math.abs(point[2]);
+          if(Math.abs(specz-abz) < VxdHitRerr)
+          {
+		   double tle = point[2]/acrad;
+		   tle+=5.0;
+		   if(tle<0.) tle=0.;
+		   int tind = (int) (tle*10.);
+		   if(tind>99) tind=99;
+		   tind+=slayer;
+		   nlayer[tind]++;
+		   ncl++;
+          }
+		  else
+		  {
+			System.out.println("EC VXD hit z position: "+abz+" while spec: "+specz);  
+		  }  
+		 } 
+	    } 

 		for (Enumeration hits = trkhits.getHits();(hits.hasMoreElements() && ncl<MaxPoints);) 
            {
 			TrackerHit hit = (TrackerHit) hits.nextElement();
             if(!hit.isEndcap())
             {
 			 int tlayer = hit.getLayer();

-             int layer = tlayer+VXDLayers;

+             int layer = tlayer+VXDLayers+VXDECLayers;

 		     int slayer = layer*100;
 			 if (Hhits) histogram("Hits/layer").fill(layer);
 			 double point[] = hit.getPoint();

@@ -371,7 +467,7 @@

             if(hit.isEndcap() && includeEC)
             {
 			 int tlayer = hit.getLayer();

-             int layer = tlayer+VXDLayers+CTRLayers;

+             int layer = tlayer+VXDLayers+VXDECLayers+CTRLayers;

 			 int slayer = layer * 100;
 			 if (Hhits) histogram("Hits/layer").fill(layer);
 			 double point[] = hit.getPoint();

@@ -392,7 +488,7 @@

             }
 		}
 

-        if(dohist)

+/*        if(dohist)

 		{
 			int nnhits = 0;
 			for(int i=0; i<100; i++)

@@ -431,7 +527,7 @@

 			for(int i=ntrst+400; i<ntrst+500; i++)
 			 nnhits += nlayer[i];
 			h025.fill(nnhits); 

-		}

+		} */

 		if (printPoints && ncl>10) System.out.println(" ...");
 	        if(firstEvents) System.out.println(" Total of "+ncl+" hits.");
 

@@ -442,7 +538,7 @@

 		{   ipoi[layer] = ipoi[layer-1]+nlayer[layer-1];
 			ilon[layer] = 0;
 		}

-        tpcsti = ipoi[VXDLayers*100];

+        tpcsti = ipoi[(VXDLayers+VXDECLayers)*100];

 	    int vhind=0;
 		for (Enumeration hits = vxdhits.getHits();(hits.hasMoreElements() && nPoints<MaxPoints);) {
 			vhind++;

@@ -450,40 +546,95 @@

 			int layer = hit.getLayer();
             if(layer > 127) layer-=128;
 			int slayer = layer*100;

-            double specr=vspecs.getLayerRadius(layer);
-			double point[] = hit.getPoint();
-            double acrad=Math.sqrt(point[0]*point[0]+point[1]*point[1]);
-			double tle = point[2]/specr;
-			tle+=5.0;
-			if(tle<0.) tle=0.;
-			int tind = (int) (tle*10.);
-			if(tind>99) tind=99;
-			tind+=slayer;
-			int iptr = ipoi[tind]+ilon[tind];
-			double dpp = Math.sqrt((point[0]-prepnt[0])*(point[0]-prepnt[0]) +
-			(point[1]-prepnt[1])*(point[1]-prepnt[1])+(point[2]-prepnt[2])*(point[2]-prepnt[2]));
-			prepnt[0]=point[0];
-			prepnt[1]=point[1];
-			prepnt[2]=point[2];
-            if((Math.abs(specr-acrad) < VxdHitRerr)&&(dpp>MinVxdHitSep))
-            {
-			 x[iptr] = point[0];
-			 y[iptr] = point[1];
-			 z[iptr] = point[2];
-			 KFlag[iptr] = false;
-             tcphit[iptr] = hit;
-			 ilon[tind]++;
-			 nPoints++;
-                        }
-		}
-		System.out.println("VXD has "+nPoints+" hits");

+			if(!hit.isEndcap())
+			{
+             double specr=vspecs.getLayerRadius(layer);
+			 double point[] = hit.getPoint();
+			 if((trace) && (mcptt != null))
+			 {
+			  if(hit.getMCParticle() == mcptt)
+			  {
+			   System.out.println
+			   ("tr. MCP hit in b. VXD lr "+layer+" x= "+df.format(point[0])+
+			   " y= "+df.format(point[1])+" z= "+df.format(point[2]));
+//			   h100.fill(point[0],point[1]);
+               double hr = Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+//               h101.fill(point[2],hr); 			   
+		      }   
+			 } 
+             double acrad=Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+			 double tle = point[2]/specr;
+			 tle+=5.0;
+			 if(tle<0.) tle=0.;
+			 int tind = (int) (tle*10.);
+			 if(tind>99) tind=99;
+			 tind+=slayer;
+			 int iptr = ipoi[tind]+ilon[tind];
+			 double dpp = Math.sqrt((point[0]-prepnt[0])*(point[0]-prepnt[0]) +
+			 (point[1]-prepnt[1])*(point[1]-prepnt[1])+(point[2]-prepnt[2])*(point[2]-prepnt[2]));
+			 prepnt[0]=point[0];
+			 prepnt[1]=point[1];
+			 prepnt[2]=point[2];
+//             if((Math.abs(specr-acrad) < VxdHitRerr)&&(dpp>MinVxdHitSep))
+             if(Math.abs(specr-acrad) < VxdHitRerr)
+             {
+			  x[iptr] = point[0];
+			  y[iptr] = point[1];
+			  z[iptr] = point[2];
+			  KFlag[iptr] = false;
+              tcphit[iptr] = hit;
+			  ilon[tind]++;
+			  nPoints++;
+             }
+		    }
+			if(hit.isEndcap() && includeVXDEC)
+			{
+		     int tlayer = layer;	 
+		     layer+= VXDLayers;
+		     slayer = layer*100;
+		     double point[] = hit.getPoint();
+			 if((trace) && (mcptt != null))
+			 {
+			  if(hit.getMCParticle() == mcptt)
+			  {
+			   System.out.println("tr. MCP hit in e. VXD lr "+tlayer+" x= "+df.format(point[0])+
+			   " y= "+df.format(point[1])+" z= "+df.format(point[2]));
+//			   h100.fill(point[0],point[1]);
+               double hr = Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+//               h101.fill(point[2],hr); 			   
+		      }   
+			 } 
+             double acrad=Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+             double specz=vespecs.getLayerZ(tlayer);
+             double abz=Math.abs(point[2]);
+             if(Math.abs(specz-abz) < VxdHitRerr)
+             {
+		      double tle = point[2]/acrad;
+		      tle+=5.0;
+		      if(tle<0.) tle=0.;
+		      int tind = (int) (tle*10.);
+		      if(tind>99) tind=99;
+		      tind+=slayer;
+			  int iptr = ipoi[tind]+ilon[tind];
+			  x[iptr] = point[0];
+			  y[iptr] = point[1];
+			  z[iptr] = point[2];
+			  KFlag[iptr] = false;
+              tcphit[iptr] = hit;
+			  ilon[tind]++;
+			  nPoints++;
+	         }   
+		    }
+		   }
+//		 System.out.println("Event has "+vhind+" VXD hits");  
+//		System.out.println("Accepted "+nPoints+" VXD hits, ind of first non vxd hit is"+tpcsti);

                 int npris = 0;
 		for (Enumeration hits = trkhits.getHits();(hits.hasMoreElements() && nPoints<MaxPoints);) {
 			MutableTrackerHit hit = (MutableTrackerHit) hits.nextElement();

-                        if(!hit.isEndcap())
-                        {

+            if(!hit.isEndcap())
+            {

 			 int tlayer = hit.getLayer();

-                         int layer = tlayer+VXDLayers;

+             int layer = tlayer+VXDLayers+VXDECLayers;

 			int slayer = layer*100;
 			double point[] = hit.getPoint();
             double specr=specs.getLayerRadius(tlayer);

@@ -495,12 +646,22 @@

 			if(tind>99) tind=99;
 			tind+=slayer;
 			int iptr = ipoi[tind]+ilon[tind];

+			if((trace) && (mcptt != null))
+			{
+			 if(hit.getMCParticle() == mcptt)
+			 {
+			   System.out.println("tr. MCP hit in b. trkr lr "+tlayer+" x= "+df.format(point[0])+
+			   " y= "+df.format(point[1])+" z= "+df.format(point[2]));
+//			   h200.fill(point[0],point[1]);
+               double hr = Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+//               h201.fill(point[2],hr); 			   
+		     }   
+			} 

             if(Math.abs(specr-acrad) < TrkrHitRerr)
             {
 			 x[iptr] = point[0];
 			 y[iptr] = point[1];
 			 z[iptr] = point[2];

-//             if(layer == 9){ z[iptr] = 1000.; point[2]=1000.; hit.setPoint(point); }

 			 KFlag[iptr] = false;
              tcphit[iptr] = hit;
 			  ilon[tind]++;

@@ -510,7 +671,7 @@

             if(hit.isEndcap() && includeEC)
             {
 			 int tlayer = hit.getLayer();

-             int layer = tlayer+VXDLayers+CTRLayers;

+             int layer = tlayer+VXDLayers+VXDECLayers+CTRLayers;

 			 int slayer = layer * 100;
 			 double point[] = hit.getPoint();
              double acrad=Math.sqrt(point[0]*point[0]+point[1]*point[1]);

@@ -523,6 +684,17 @@

 			 if(tind>99) tind=99;
 			 tind+=slayer;
 			 int iptr = ipoi[tind]+ilon[tind];

+			if((trace) && (mcptt != null))
+			{
+			 if(hit.getMCParticle() == mcptt)
+			 {
+			   System.out.println("tr. MCP hit in ec trkr lr "+tlayer+" x= "+df.format(point[0])+
+			   " y= "+df.format(point[1])+" z= "+df.format(point[2]));
+//			   h200.fill(point[0],point[1]);
+               double hr = Math.sqrt(point[0]*point[0]+point[1]*point[1]);
+//               h201.fill(point[2],hr); 			   
+		     }   
+			} 

              if(Math.abs(specz-abz) < EcHitZerr)
              {
 			  x[iptr] = point[0];

@@ -535,11 +707,12 @@

              }
             }
 		}

-		System.out.println("Total number of hits in tracking system "+nPoints);
-		for (int level=MinLevel; level<=MaxLevel; level++) {
-			if (first) System.out.println(" Find tracks at level "+level);
-			finder.findTracks(level,nPoints,NLayers);
-			if (firstEvents)

+//		System.out.println("Total number of hits in tracking system "+nPoints);
+		for (int level=MinLevel; level<=MaxLevel; level++) 
+		{
+		 if (first) System.out.println(" Find tracks at level "+level);
+		 finder.findTracks(level,nPoints,NLayers);
+		 if (firstEvents)

 				System.out.println("  Number of tracks at this level = "+finder.getNTracks());
 		}
 		// Get number of tracks found

@@ -560,7 +733,7 @@

           for(int j=1; j<=ntp; j++)
              if(itra[j]<tpcsti) nvxpnt[i]++;
 		  VD.assignHits(i);				 

-			if (Hstd) histogram("NPoints/track").fill(getNPoints(i));

+		  if (Hstd) histogram("NPoints/track").fill(getNPoints(i));

 			if (Htrks) {
 				histogram("d0").fill(finder.get_d0(i));
 				histogram("phi0").fill(finder.get_phi0(i));

lcd/hep/lcd/recon/tracking/ccd

diff -u -r1.1 -r1.2
--- VXDRecoVD.java	11 Mar 2005 19:41:00 -0000	1.1
+++ VXDRecoVD.java	10 Aug 2005 20:28:17 -0000	1.2
@@ -28,8 +28,11 @@

 			   implements VertexDetector
 {
     Specifications specs;

+    Specifications especs;

     public boolean trace = false;
 	private boolean initialized = false;

+    int Nvl = 5;
+	int Nevl = 4;

     int NLayers=5;
     int MaxPoints = 10000;
     int MaxTracks = 2000;

@@ -52,6 +55,7 @@

     {
 	 this.tracker = trckr;
      specs = tracker.getVxdSpecifications();

+	 especs = tracker.getVxdECSpecifications();

     }
 
     public Specifications getSpecifications()

@@ -141,12 +145,15 @@

 
     public void startEvent(LCDEvent event)
     {

-      specs.update(event);
-	  if(specs.getNLayers() != NLayers) initialized = false;

+	  Nvl = specs.getNLayers();
+	  Nevl = especs.getNLayers();
+	  int ntvl = Nvl;
+	  if(tracker.includeVXDEC) ntvl = Nvl+Nevl;
+	  if(ntvl != NLayers) initialized = false;

 	  if(!initialized)
 	  {
 	   finder = (VXDPatFind) (tracker.getFinder());	  

-	   NLayers = specs.getNLayers();	  

+	   NLayers = ntvl;	  

 	   MaxTrkPts = NLayers;	  
 	   ipnt = new int[MaxTracks][MaxTrkPts+1];
 	   MaxPoints = MaxTracks * MaxTrkPts;

@@ -193,6 +200,11 @@

 		int layer = hit.getLayer();
 		if(layer < MaxTrkPts) order[layer]=i;
 	   }   

+	   if((hit.getSystem() == hit.SYSTEM_VXD) && (hit.isEndcap()))
+	   {
+		int layer = hit.getLayer()+ Nvl;
+		if(layer < MaxTrkPts) order[layer]=i;
+	   }   

 	  }  
 	 }
 	 for(int i=0; i<MaxTrkPts; i++)