////////////////////////////////////////////
      ////////// FINAL correctionn ////////////////
      ///////////////////////////////////////////
      
      TCanvas *cxz = new TCanvas("cxz","hit pattern (slope corrected)",1010,10,510,510);
      gr ->Draw("AC");
      gs ->Draw("same");
      gt ->Draw("same");
      //cell boundary
      Float64_t xr[2], yr[2];
      for (Int_t k=0;k<12;k++) {
	float64 ang=(30.*k)/57.296;
	xr[0] = RII*cos(ang);
	yr[0] = RII*sin(ang);
	xr[1] = RO*cos(ang);
	yr[1] = RO*sin(ang);
	TGraph *gg = new TGraph(2,xr,yr);
	gg ->SetLineStyle(3);
	gg ->SetLineWidth(2);
	gg ->SetLineColor(2);
	gg ->Draw("same");
      }
      for (Int_t k=0;k<12;k++) {
	float64 ang=(15.0+30.*k)/57.296;
	xr[0] = RII*cos(ang);
	yr[0] = RII*sin(ang);
	xr[1] = RO*cos(ang);
	yr[1] = RO*sin(ang);
	TGraph *gh = new TGraph(2,xr,yr);
	gh ->Draw("same");
      }
      cout <<" "<< endl;
      
      Float64_t fxh[40], fyh[40];
      for(k=1 ; k<40;k++){
	fxh[k]=-10000;
	fyh[k]=-10000;}
      int hit2channel[10];
      for (int w=0; w<10; w++){
	hit2channel[w]=100;}
      
      int iterend=10;  //two more iterations
      float64 rmsMax=10000;
      for(int iteration=3; iteration <(iterend+1) ;iteration++){
	cout<<"Final iteration"<<iteration<<"  "<<newslp<<" SI "<<newint<< endl;
	int khit=0;
	float64 xbottom=-2.5/newslp-newint/newslp;
	float64 xtop   = 2.5/newslp-newint/newslp;
	float64 xb=xbottom;
	float64 xe=xtop;
	if(xtop <xbottom) {
	  xb=xtop;
	  xe=xbottom;     }
	float64 xstep=(xe-xb)/1000.;
	//cout <<" x12: "<<xb<<" "<<xe<<" "<<xstep<< endl;
	
	for(int jj=0;jj<16;jj++){
	  if(hit[jj] == 0) continue;
	  float64 time=Ltime[1][jj]-t0;
	  float64 phimin=30.0*(jj-1)/rad;
	  float64 phimax=phimin+30.0/rad;
	  
	  //step through the cell following track
	  float64 rnow  =10000.0;
	  float64 phinow=10000.0;
	  float64 tnow  =10000.0;
	  int jrnow=-100;
	  int jphinow=-100;
	  //Follow the track
	  float64 Tfromtable=10000.0;
	  float64 EX=10000.0;
	  float64 EY=10000.0;
	  
	  for(float64 xj=xb; xj< xe ;xj=xj+xstep){
	    float64 yj=newslp*xj+newint;
	    float64 rj=sqrt(xj*xj+yj*yj);
	    float64 phij=aatan(yj,xj);
	    //cout << jj<<" in  "<<xj<<" "<<yj<<" "<<rj<<" "<<phij<< endl;
	    if(phij < phimin || phij > phimax) continue;
	    if(rj  < 0.15    || rj  > 2.75   ) continue;
	    int jr=rj/0.01;
	    int jp=(phij*rad-30*(jj-1))/0.5;
	    //cout << jj<<" out "<<xj<<" "<<yj<<" "<<rj<<" "<<phij<<" ("<<jr<<" "<<jp<<")"<< endl;
	    float64 Tfromtable=Tstore[jr][jp];
	    //cout <<jr<<" "<<jp<<" "<<Tfromtable<< endl;
	    if(Tfromtable >500.0) continue;
	    if(Tfromtable > tnow) continue;
	    tnow=Tfromtable;
	    phinow=phij;
	    rnow=rj;
	    jrnow=jr;
	    jphinow=jp;
	    //EX=EXstore[jr][jp];
	    //EY=EYstore[jr][jp];
	  }
	  if(rnow > 2.65) continue;
	  float64 xnow=rnow*cos(phinow);
	  float64 ynow=rnow*sin(phinow);
	  float64 angle=wireAngle[jj];
	  float64 xwire=RI*cos(angle);
	  float64 ywire=RI*sin(angle);
	  float64 xprime=xnow-xwire;
	  float64 yprime=ynow-ywire;
	  float64 rprime=sqrt(xprime*xprime+yprime*yprime);
	  float64 rnew=rprime+0.005*(time-tnow); //correct for the time difference
	  if(rnew > 1.25) rnew=1.25;
	  float64 xnew= xwire+rnew*xprime/rprime;   
	  float64 ynew= ywire+rnew*yprime/rprime;   
	  //cout <<jj<<" "<<khit<<" "<<xnew<<" :rphimin1: "<<ynew<<" ("<<tnow<<" :t: "<<time<<")"<< endl;
	  float64 rot=angle-wireAngle[1];
	  EX=1000.0; 
	  EY=268.0; 
	  if(rnow > RI) { 
	    EX=-1000.0; 
	    EY=-268.0;  } 
	  float64 EXr= EX*cos(rot)-EY*sin(rot);
	  float64 EYr= EX*sin(rot)+EY*cos(rot);
	  float64 ET=sqrt(EXr*EXr+EYr*EYr);
	  float64 xxnew=xwire+xprime-0.005*(time-tnow)*EXr/ET;
	  float64 yynew=ywire+yprime-0.005*(time-tnow)*EYr/ET;
	  float64 rrnew=sqrt(xxnew*xxnew+yynew*yynew);
	  if(rrnew < 0.23 || rrnew > 2.52) continue;
	  //cout <<jj<<" "<<khit<<" "<<xxnew<<" :rphimin2: "<<yynew<<" ("<<tnow<<" :t: "<<time<<")"<< endl;
	  //cout <<"  R: "<<rot<<" "<<EX<<" "<<EY<<" "<<EXr<<" "<<EYr<< endl;
	  fxh[khit]=xxnew;
	  fyh[khit]=yynew;
	  hit2channel[khit]=jj;
	  khit=khit+1;
	}
	//TGraphErrors *gfinal = new TGraphErrors(khit,fyh,fxh,yerror,xerror);
	TGraph *gfinal = new TGraph(khit,fyh,fxh);
	TGraph *gfdata = new TGraph(khit,fxh,fyh);  //this is for plot
	gfdata->SetMarkerColor(4);
	gfdata->SetMarkerStyle(21);   
	gfdata->SetMarkerSize(2);   
	pol1 ->SetParameter(0,-newint/newslp);
	pol1 ->SetParameter(1,1.0/newslp);
	gfinal ->Fit("pol1","Q");
	float64 chi2 = pol1 ->GetChisquare();
	float64 newslp = pol1 ->GetParameter(1);
	float64 newint = pol1 ->GetParameter(0);
	newint=-newint/newslp;
	newslp=1.0/newslp;
	
	float64 rmsSaved[10];
	int jjdelete=0;
	int Nlayer=khit;
	float64 fxq[10],fyq[10];
	float64 sumdely=0.0;
	float64 sum2dely=0.0;
	for(int w; w<10; w++){
	  rmsSaved[w]=-1000.;}
	for (int jk=0; jk<Nlayer; jk++){
	  fyq[jk]=fyh[jk];
	  fxq[jk]=1.0/newslp*fyq[jk]-newint/newslp;
	  float64 dely=10000.*(fxq[jk]-fxh[jk]);
	  sumdely =sumdely+dely;
	  sum2dely=sum2dely+dely*dely;
	  rmsSaved[jk]=fabs(dely);}
	
	float64 rms=sqrt(sum2dely/float(Nlayer));
	cout<<"stat : N= "<<Nlayer<<",  RMS= "<<rms<<" microns"<< endl;    
	// if(rms < 800.0 && (iteration == iterend || rms < 100. || (rms <rmsMax && (rmsMax-rms) <10.0))) {
	// if(rms < 1000.0 && (iteration == iterend || rms < 100.)) {
	if(rms < 1000.0 && iteration == iterend) {
	  gfdata ->Draw("same*");  
	  float64 fxa[600],fya[600];
	  for (int jk=0; jk<600; jk++){
	    fxa[jk]=-3.0+0.01*jk;
	    fya[jk]=newslp*fxa[jk]+newint;}
	  TGraph *gsw = new TGraph(600,fxa,fya);
	  gsw->SetLineColor(2);
	  gsw->SetLineWidth(2);
	  gsw ->Draw("same");
	  sprintf(buffer,"rms:%d um",rms);
	  TLegend *lp1 = new TLegend(0.65,0.80,1.1,0.9);
	  lp1->SetFillStyle(0);
	  lp1 ->SetBorderSize(0);
	  lp1->SetHeader(buffer);
	  lp1 ->Draw();
	  break;      }
	rmsMax=rms;
	
	if(iteration == 5) {
	  float64 rmsdelete=0.0;
	  int jjdelete=0;
	  for(int w=0; w<10;w++){
	    if(rmsSaved[w] < 0.0) continue;
	    if(rmsSaved[w] < rmsdelete) continue;
	    rmsdelete=rmsSaved[w];
	    jjdelete=w;         }
	  if(rmsdelete < 500.0) continue;
	  cout <<"*hit deleted: "<<rmsdelete<<endl;
          //delete the hit and refit    
          hit[hit2channel[jjdelete]]=0;
	}
	
      } //end final iteration
      
      cxz ->Update();
    ENDPLOT2:
      cxx ->Update();
      cxy ->Update();
      
      cout <<"DONE2! Event number = "<<eventNow<<" @ "<<i<< endl;
      int more;
      cout <<" "<<endl;
      cout <<" "<<i<<" "<<vsize << endl;
      cout <<"  More?? " << endl;
      cin >> more;  // pause for the next event.
      
      if(more == -1) goto REPEAT;      
      for(int jj=0;jj<40;jj++){
	hit[jj]=0;
	for(int jk=0;jk<10;jk++){
	  Ltime[jk][jj]=0;
	  Ttime[jk][jj]=0;
	  Width[jk][jj]=0;
	}}
      
      if(more < 0) {
	eventPlot = abs(more);
	goto BEGIN;  }