#include <fcntl.h>
#include <sys/ioctl.h>
#include <stdio.h>
#include <math.h>

float chifit_ts(float tql,float tq,float tqh);
float chifit_gatti(float xa,float ts,float *err,float *d,int *u);

extern float xalct;
extern int alctgood,alctbad; 

void cat_trk_find(int ievt,int brd);

extern unsigned short int data[16][16][6][5];
int use_pt[6];
float dpt[6][9];
float x_pt[6],y_pt[6],dy_pt[6],qchi_pt[6],gat_pt[6];
float ts[6],max[6];
float offset[6]={-0.002,-0.493,-0.005,-0.496,0.001,-0.499};

void cat_trk_find(int ievt,int brd)
{
int layer[6]={2,0,4,5,3,1};
int maxval,maxch,maxtim;
float cgatti;
int i,j,k,l,m,n;
float d[9];
int shit;
int u[3];
float err_xx[3];

/* crude track finding - just take maximum in each plane */
  for(i=0;i<6;i++){
    maxval=0;
    for(j=0;j<16;j++){
      for(k=0;k<16;k++){
	if(data[k][j][i][brd]-data[0][j][i][brd]>maxval){
          maxval=data[k][j][i][brd]-data[0][j][i][brd];
          maxch=j;
          maxtim=k;
        }
      }
    }     
    // printf(" plane %d %d maxch %d maxtim %d \n",layer[i],i,maxch,maxtim);
    if(maxch<2||maxch>13)goto Loopr;
    n=0;
    for(l=-1;l<2;l++){
      for(m=-1;m<2;m++){
      dpt[layer[i]][n]= data[maxtim+l][maxch+m][i][brd]-(data[0][maxch+m][i][brd]+data[1][maxch+m][i][brd])/2.; 
      n=n+1;
      }
    }
    
/* now determine position of points */
 
    /* fit for start time */

   ts[layer[i]]=chifit_ts(dpt[layer[i]][1],dpt[layer[i]][4],dpt[layer[i]][7]);     

   /* fit for position */
   
   for(l=0;l<9;l++)d[l]=dpt[layer[i]][l];  
   u[0]=1;u[1]=1;u[2]=1;
   cgatti=chifit_gatti(xalct,ts[layer[i]],err_xx,d,u);

   use_pt[layer[i]]=1;
   x_pt[layer[i]]=layer[i];
   y_pt[layer[i]]=maxch+cgatti+offset[layer[i]];
   dy_pt[layer[i]]=err_xx[0];
   qchi_pt[layer[i]]=err_xx[2];
   gat_pt[layer[i]]=cgatti;
   // printf(" i %d %f %f %f %f \n",i,x_pt[layer[i]],y_pt[layer[i]],dy_pt[layer[i]],qchi_pt[layer[i]]);
   
  Loopr:        

  }

}