public void FitTracks(double _r[], double _phi[], double _z[],
						double _dr[], double _dphi[], double _dz[],
						int _layern[], int nhits, Appearance trackAppearance) {

		final int    MAXLAYER = 40;                  // no of layers in the tracker

		final double CHISQMAX = 200.;                // Chisq minimum beyond which a global fit is bad 
		final double CHISQSEED = 100.;               // Chisq minimum for the three seed hit helix (should be more stringent)
		final double CHISQMAXHIT = 20.;              // Chisq maximum for an individual hit                          
		final double MAXDPHI = 0.2*Math.PI;          // Maximum difference in phi between candidate seeds and hits on tracks
													 // Note that phi runs from -PI to +PI

		final double CUTR = 20.;                     // Used to check intercept at R axis for two seed hits
		final double STRIPWIDTH = 10.;               // Used for simple cut when extrapolating track to hit
		final double MINRTR = 84.;                   // minimum acceptable radius of helix fit
													 // Note that 100.Pt(Gev)/(0.3.B(Tesla)) = Rtr(cm.) thus 1Gev ~= 84 cm.
		final int    MINGOODHITS = 3;                // Minimum number of well-measured hits on a track
		final int    MINHITS = MINGOODHITS+3;        // Minimum number of hits on a track

		int lhit [][] = new int[nhits][MAXLAYER];    // index into hit list ordered by layer
		int ninlayer [] = new int[MAXLAYER];         // Number of hits in each layer
		boolean used [] = new boolean[nhits];        // flag to indicate hit is used in a track fit
		double rhighlayer [] = new double[MAXLAYER]; // maximum R coordinate found on each layer


		long starttime, endtime;

		for (int i=0;i<MAXLAYER;i++) {
			ninlayer[i] = 0;
			rhighlayer[i] = 0.;
		}
		for (int i=0;i<nhits;i++) {
			used[i] = false;
			int ilay = _layern[i];
			int ninl = ninlayer[ilay];
			lhit[ninl][ilay] = i;
			ninlayer[ilay]++;
			if(_r[i] > rhighlayer[ilay]) rhighlayer[ilay] = _r[i];
		}

		// Find trios of hits to seed a track

		double [] r   = new double[3];
		double [] phi = new double[3];
		double [] z   = new double[3];
		double [] dr  = new double[3];
		double [] dphi= new double[3];
		double [] dz  = new double[3];

		int nUsedSoFar = 0; // tally of number of hits accounted for in track fits
		int ntracks = 0;
		starttime = System.currentTimeMillis(); // Mark the current system time

OUTER_LAYER:
		for (int layer1=MAXLAYER-1;layer1>1;layer1--) { // OUTER LAYER
OUTER_LAYER_HITS:
			for (int i1=0;i1<ninlayer[layer1];i1++) {
				int ih1 = lhit[i1][layer1];
				if(used[ih1]) continue; // select only unused hits
				dr[0]   = _dr[ih1];
				dphi[0] = _dphi[ih1];
				dz[0]   = _dz[ih1];
				if(dr[0] > 1. || dz[0] > 1.) continue; // select only stereo or pixel hits
				r[0]    = _r[ih1];
				phi[0]  = _phi[ih1];
				z[0]    = _z[ih1];

				//System.out.println(layer1);

MIDDLE_LAYER:
				for (int layer2=layer1-1;layer2>0;layer2--) { 
					for (int i2=0;i2<ninlayer[layer2];i2++) {
						int ih2 = lhit[i2][layer2];
						if(used[ih2]) continue; // select only unused hits
						dr[1]   = _dr[ih2];
						dphi[1] = _dphi[ih2];
						dz[1]   = _dz[ih2];
						if(dr[1] > 1. || dz[1] > 1.) continue; // select only stereo or pixel hits
						r[1]    = _r[ih2];
						phi[1]  = _phi[ih2];
						z[1]    = _z[ih2];

						double DPHI = Math.abs(phi[0] - phi[1]);
						if(DPHI > Math.PI) DPHI = 2.*Math.PI - DPHI;
						if(DPHI > MAXDPHI) continue;

						// Test that this pair and the centre of the detector make a valid straight line fit
						dr[2]   = 1.;
						dphi[2] = MAXDPHI; 
						dz[2]   = 1.;
						r[2]    = 0.;
						phi[2]  = 0.5*(phi[1] + phi[0]); // Very abitrary but we're only making a rough check
						z[2]    = 0.;

						// Fit a straight line in the RZ plane
						double a=0.;
						double b=0.;
						double chisq = FitLine(phi,dphi,r,dr,z,dz,a,b);
						//System.out.println((float)chisq+" "+a);
						if(chisq > 3.*CHISQMAXHIT) continue;
						System.out.print("|");
						double [] parameters = new double[7];
INNER_LAYER:
						for (int layer3=0;layer3<layer2;layer3++) { // INNER LAYER
INNER_LAYER_HITS:
							for (int i3=0;i3<ninlayer[layer3];i3++) {
								int ih3 = lhit[i3][layer3];
								if(used[ih3]) continue; // select only unused hits
								dr[2]   = _dr[ih3];
								dphi[2] = _dphi[ih3];
								dz[2]   = _dz[ih3];
								if(dr[2] > 1. || dz[2] > 1.) continue; // select only stereo or pixel hits
								r[2]    = _r[ih3];
								phi[2]  = _phi[ih3];
								z[2]    = _z[ih3];

								DPHI = Math.abs(phi[1] - phi[2]);
								if(DPHI > Math.PI) DPHI = 2.*Math.PI - DPHI;
								if(DPHI > MAXDPHI) continue;

								// Fit a helix, asking for a min. Rtr of MINRTR cm.
								chisq = FitHelix(phi,dphi,r,dr,z,dz,MINRTR,parameters);
								if(chisq > CHISQSEED) continue;

								int nHits = 0;
								int nGoodHits = 0;
								int trackHits [] = new int[MAXLAYER];
								boolean trackGoodHits [] = new boolean[MAXLAYER];
								double trackChisq = 0.;

								// Candidate track. Move outwards from innermost layer marking closest hits

								for (int layer4=0;layer4<MAXLAYER;layer4++) {
									if(ninlayer[layer4] <= 0) continue;
									// Make sure it makes sense to check this layer
									double rTrack = Math.sqrt(parameters[1]*parameters[1]+parameters[2]*parameters[2]);
									if (rTrack > rhighlayer[layer4]) continue;
									double zTrack = parameters[0];
									double phiTrack = Math.atan2(parameters[2],parameters[1]);
									int best = -1;
									double chimin = 99999.;
									double dchisq;
									double bestpar [] = new double[7];
									double oldpar [] = new double[7];
									for (int i=0;i<7;i++) oldpar[i] = parameters[i];
									for (int i4=0;i4<ninlayer[layer4];i4++) { 
										int ih4 = lhit[i4][layer4];
										if(used[ih4]) continue;
										// Extrapolate the track to the hit's Z or R position 
										double Z    = _z[ih4];
										double R    = _r[ih4];
										double Phi  = _phi[ih4];
										// Check that the difference in Phi is not wild
										DPHI = Math.abs(phiTrack - Phi);
										if (DPHI > Math.PI) DPHI = 2.*Math.PI - DPHI;
										if(DPHI > MAXDPHI) continue;
										// Check that the Z,R of the hit is on the correct side
										if(zTrack > 0. && Z < 0.) continue;
										if(zTrack < 0. && Z > 0.) continue;
										if(rTrack-R > STRIPWIDTH) continue; // Safe for strip width
										if(zTrack > 0. && zTrack-Z > STRIPWIDTH) continue;
										if(zTrack < 0. && Z-zTrack > STRIPWIDTH) continue;
										double dZ   = _dz[ih4];
										double dPhi = _dphi[ih4];
										double dR   = _dr[ih4];
										for(int i=0;i<7;i++) parameters[i] = oldpar[i];
										if(dR > dZ) {
											dchisq = ExtrapolateToPlane(Phi,dPhi,R,dR,Z,dZ,parameters);
										} else {
											dchisq = ExtrapolateToCylinder(Phi,dPhi,R,dR,Z,dZ,parameters);
										}
										if(dchisq < chimin) {
											best = ih4;
											chimin = dchisq;
											for (int i=0;i<7;i++) bestpar[i] = parameters[i];
										}
									}
									if(best == -1 || chimin >= CHISQMAXHIT) {
										for (int i=0;i<7;i++) parameters[i] = oldpar[i];
										continue;
									} else {
										trackChisq += chimin;
										//System.out.println("layer dchi="+chimin+" current track chi="+trackChisq);
										if(trackChisq > CHISQMAX) {
											// track must be rejected
											// continue with the next hit in the three seed hit sequence
											continue INNER_LAYER_HITS;
										}
										for(int i=0;i<7;i++) parameters[i] = bestpar[i];
										trackHits[nHits] = best;
										if(_dr[best] < 1. && _dz[best] < 1.) {
											trackGoodHits[nHits] = true;
											nGoodHits++;
										} else {
											trackGoodHits[nHits] = false;
										}
										nHits++;
									}
								}
								// Finished looping over all layers for this candidate track
								//System.out.println("Finished layer loop :"+nHits+" "+nGoodHits);
								if(trackChisq > CHISQMAX || nGoodHits < MINGOODHITS || nHits < MINHITS) {
									continue INNER_LAYER_HITS; // return for next inner layer hit in seed trio
								} else {
									// Good track
									// Fit again using up to the first three good hits on the track
									int in = 2;
									int goodhitsleft = nGoodHits;
									for(int i=0;i<nHits;i++) {
										if(nHits-i > in+1) {
											if(goodhitsleft > 0 && !trackGoodHits[i]) continue;
										}
										if(trackGoodHits[i]) goodhitsleft--;
										//if(trackGoodHits[i]) System.out.print("-good-");
										//if(!trackGoodHits[i]) System.out.print("-bad-");
										int ip = trackHits[i];
										r[in] = _r[ip];
										phi[in] = _phi[ip];
										z[in] = _z[ip];
										dr[in] = _dr[ip];
										dphi[in] = _dphi[ip];
										dz[in] = _dz[ip];
										if(in <= 0) break;
										in--;
									}
									if(in != 0) continue INNER_LAYER_HITS; // Sanity

									// Fit a helix, with a min. Rtr of MINRTR cm.
									chisq = FitHelix(phi,dphi,r,dr,z,dz,MINRTR,parameters);
									if(chisq < CHISQMAX) {
										System.out.println("\nBest chisq="+trackChisq);
										System.out.println("Pt="+0.01*0.3*4.0/parameters[5]);
										System.out.println(nHits+" points on track ("+nGoodHits+" well measured)");
										System.out.println("Refit first well measured 3 points: chisq="+chisq);
										System.out.println(parameters[0]+" "+parameters[1]+" "+parameters[2]+" "+parameters[3]+" "+parameters[4]+" "+parameters[5]);
										ntracks++;
										Track3D Track3DObj = new Track3D(parameters,trackAppearance);
										Shape3D track = Track3DObj.get();
										//track.setCapability(Shape3D.ALLOW_APPEARANCE_READ);
										//track.setCapability(Shape3D.ALLOW_APPEARANCE_WRITE);
										if(track != null) {
											eventBG.addChild(track);
										}
										// mark points as used
										for(int i=0;i<nHits;i++) {
											nUsedSoFar++;
											used[trackHits[i]] = true;
										}
										//if(ntracks >= 4) break OUTER_LAYER;
										if(nhits-nUsedSoFar < 3) {	
											break OUTER_LAYER; // Fitting has finished since there are no hits left to use
										} else {
											continue OUTER_LAYER_HITS; // Go for a completely fresh new 3 hit seed
										}
									}

								} // INNER_LAYER
							}
						}
					}
				}
			}
		}
		endtime = System.currentTimeMillis();
		int dtime = (int) (endtime-starttime);

		System.out.println("\nFitting complete for "+nhits+" hits after "+dtime+" milliseconds");
		System.out.println("There were "+ntracks+" tracks found, and "+nUsedSoFar+" hits used in the fits");
	}



	public double FitHelix(double phi[], double dphi[], 
		                 double r[], double dr[],
						 double z[], double dz[],
						 double minRtr,
						 double parameters[]) {
		// This function takes as input 3 points (r,phi,z) and a minimum acceptable
		// radius of curvature (minRtr), and attempts to fit a helix to the points.
		// The function returns the chisq of the fit as its value, and the
		// fitted parameters (if any) in parameters[], where the values are:
		// [0] = Z
		// [1] = X  coordinates of the track at point 3
		// [2] = Y
		// [3] = Theta, the dip angle, where Tan(theta) = dS/dZ i.e. the track length per unit change in Z
		// [4] = Phi, the angle in the XY plane that point 3 makes with the centre of the helix circle
		// [5] = 1/Rtr, where Rtr is the radius of the helix circle in the XY plane
		// [6] = Sin(beta) where beta is the difference between Phi and Atan(Y/X)

		double chi_squared=0.;

		double [] x = new double[3];
		double [] y = new double[3];
		// Calculate positions of the points in the XY plane
		for(int i=0;i<3;i++) {
			x[i] = r[i]*Math.cos(phi[i]);
			y[i] = r[i]*Math.sin(phi[i]);
		}
		// Calculate distances between the points in the XY plane
		double dy21 = y[1]-y[0];
		double dx21 = x[1]-x[0];
		double d21sq = dx21*dx21 + dy21*dy21;

		double dy23 = y[1]-y[2];
		double dx23 = x[1]-x[2];
		double d23sq = dx23*dx23 + dy23*dy23;

		// Calculate (signed) curvature of the line joining the points
		double temp = dx21*dy23 - dx23*dy21;
		double tempsign = 1.;
		if(temp < 0.) tempsign = -1.;

		if(Math.abs(temp) < 1.0e-10) {
			//System.out.println("Very large Rtr ... ");
			temp = tempsign*0.00001;
		} 

		double sdet = 0.5/temp;
		// Calculate the position of the centre of the helix in the XY plane
		double xc = x[1] + sdet*(d23sq*dy21 - d21sq*dy23);
		double yc = y[1] + sdet*(d21sq*dx23 - d23sq*dx21);
		double xd = x[1] - xc;
		double yd = y[1] - yc;

		double Rtr = tempsign*Math.sqrt( xd*xd + yd*yd );

		if(Math.abs(Rtr) < minRtr) return 99999.;

		double Rinv = 1./Rtr;

		// angles at point 3 and curvature distance (distance along helix) between 3-2
		double sina = (x[2]-xc)*Rinv;
		double cosa = (yc-y[2])*Rinv;

		double sinb = (sina*x[2] - cosa*y[2]) / r[2];
		if(Math.abs(sinb) >= 1.) {
			System.out.println("sinb >= 1 ... error ");
			return (99999.);
		}

		double alph3 = Math.atan2(sina,cosa);
		double alph2 = Math.atan2( (x[1]-xc)*Rinv, (yc-y[1])*Rinv );
		//double alph1 = Math.atan2( (x[0]-xc)*Rinv, (yc-y[0])*Rinv );

		double curv3 = Rtr*(alph2-alph3);

		double theta = Math.atan2(curv3,z[1]-z[2]);
		if(theta < 0.) theta += Math.PI;


		parameters[0] = z[2];
		parameters[1] = x[2];
		parameters[2] = y[2];
		parameters[3] = theta;
		parameters[4] = phi[2];
		parameters[5] = Rinv;
		parameters[6] = sinb;

		// Calculate the Chi squared of the fit by extrapolating the helix
		// from point 3 to point 2 and then point 1 and summing the chisq contributions
		double partest [] = new double[7];
		for(int i=0;i<7;i++) partest[i] = parameters[i];
		double sumchisq = 0.;
		for(int i=1;i>=0;i--) {
			if(dr[i] > dz[i]) {
				sumchisq += ExtrapolateToPlane(phi[i],dphi[i],r[i],dr[i],z[i],dz[i],partest);
			} else {
				sumchisq += ExtrapolateToCylinder(phi[i],dphi[i],r[i],dr[i],z[i],dz[i],partest);
			}
		}
		return sumchisq;
	}


	public double ExtrapolateToPlane(double phi, double dphi, 
									double r, double dr, 
									double z, double dz,
									double parameters[]) {

		double Z = parameters[0];
		double X = parameters[1];
		double Y = parameters[2];
		double tanl = Math.tan(parameters[3]);
		double PhiT = parameters[4];
		double Rinv = parameters[5];
		double sinb = parameters[6];

		double Phi = Math.atan2(Y,X);
		double R = Math.sqrt(X*X+Y*Y);

		double DZ = z-Z;
		double DS = DZ*tanl;
		double RTRK = 1./Rinv;
		double COSF0 = Math.cos(PhiT);
		double SINF0 = Math.sin(PhiT);
		double XC = X - RTRK*SINF0;
		double YC = Y + RTRK*COSF0;
		double DPHI = Rinv*DS;
		double PHI1 = PhiT + DPHI;
		if(PHI1 < 0.) PHI1 += 2.*Math.PI;
		if(PHI1 > 2.*Math.PI) PHI1 -= 2.*Math.PI;
		double COSF1 = Math.cos(PHI1);
		double SINF1 = Math.sin(PHI1);
		double X1 = XC + RTRK*SINF1;
		double Y1 = YC - RTRK*COSF1;
		double phif = Math.atan2(Y1,X1);
		double R1 = Math.sqrt(X1*X1 + Y1*Y1);

		parameters[0] = z;
		parameters[1] = X1;
		parameters[2] = Y1;
		parameters[4] = PHI1;

		//System.out.println("From track at "+R+" "+Phi+" "+Z);
		//System.out.println("Goto r,phi,z  "+r+" "+phi+" "+z);
		//System.out.println("Result        "+R1+" "+phif+" "+z);

		double dchisq = (phi-phif)*(phi-phif)/(dphi*dphi);
		dchisq += (r-R1)*(r-R1)/(dr*dr);
		//System.out.println("Z plane Chisq "+dchisq);

		return dchisq;
	}

	public double ExtrapolateToCylinder(double phi, double dphi, 
									double r, double dr, 
									double z, double dz,
									double parameters[]) {

		double Z = parameters[0];
		double X = parameters[1];
		double Y = parameters[2];
		double tanl = Math.tan(parameters[3]);
		double PhiT = parameters[4];
		double Rinv = parameters[5];
		double sinb = parameters[6];

		double R = Math.sqrt(X*X+Y*Y);

		double beta = Math.asin(sinb);
		double cosb = Math.cos(beta);
		double cotth = 1./tanl;
		double rtrk = 1./Rinv;

		double cosf0 = Math.cos(PhiT);
		double sinf0 = Math.sin(PhiT);
		double xc = X - rtrk*sinf0;
		double yc = Y + rtrk*cosf0;

		double rc2 = xc*xc + yc*yc;

		double rrr = (r*r - rtrk*rtrk - rc2) / (2.*rtrk);
		double delt = rc2 - rrr*rrr;

		if(delt <= 0.) {
			//System.out.println("Error in extrapolation");
			return 999999.;
		}

		delt = Math.sqrt(delt);

		double sinf = (xc*rrr + yc*delt)/rc2;
		double cosf = (xc*delt - yc*rrr)/rc2;

		double X1 = xc + rtrk*sinf;
		double Y1 = yc - rtrk*cosf;

		double sinbf = (sinf*X1 - cosf*Y1)/r;
		if(sinbf > 0.9999) {
			System.out.println("Beta too large at intersection");
			return 999999.;
		}

		double PHI1 = Math.atan2(sinf,cosf); // final phi around helix circle
		double dPHI = PHI1 - PhiT;           // change in phi around helix circle

		double zf = Z + cotth*rtrk*dPHI;     // resulting position in z
		
		double rf = Math.sqrt(X1*X1 + Y1*Y1);// final r in global coords
		double phif = Math.atan2(Y1,X1);     // final phi in global coords

		parameters[0] = zf;
		parameters[1] = X1;
		parameters[2] = Y1;
		parameters[4] = PHI1;
		parameters[6] = sinbf;

		//System.out.println("From track at "+R+" "+Math.atan2(Y,X)+" "+Z);
		//System.out.println("Goto r,phi,z  "+r+" "+phi+" "+z);
		//System.out.println("Result        "+rf+" "+phif+" "+zf);

		double dchisq = (phi-phif)*(phi-phif)/(dphi*dphi);
		dchisq += (z-zf)*(z-zf)/(dz*dz);
		//System.out.println("R surface Chisq "+dchisq);

		return dchisq;
	}