#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define Hlength 1326
#define Blength 658
#define Nstations 5

#define Level 100 /* 10 sigma: 1st few points fall into a single bit */

  unsigned char *buffer;
  unsigned char *Value;


int main(int argc, char **argv)
{
  FILE *fp;

  int i, j, ind;
  char *Fname;
  long lSize;

  size_t result;

  int XYstat;


/* Header variables */

  if (argc!=2) {fputs ("enter file name\n", stderr); exit(0);}

  Fname = argv[1];

  fp = fopen (Fname, "rb");
  if (fp==NULL) {fputs ("File error", stderr); exit (1);}

  /* obtain file size */

  fseek (fp , 0 , SEEK_END);
  lSize = ftell (fp);
  rewind (fp);

  /* allocate memory to contain one event */

  buffer = (unsigned char*) malloc(sizeof(unsigned char)*Hlength);

  if (buffer==NULL) {fputs ("Memory error", stderr); exit (2);}


/* copy part of the header into the buffer */

  result = fread (buffer, 1, 46, fp);
  if (result != 46) {fputs ("Reading error", stderr); exit (3);}

  int HDR[46];

  Value = buffer;

  for (i=0; i<46; i++) HDR[i] = (int)*Value++;

  int Stat, Run, RunNum[2];
  int TRG[4], Nevt;

  int COI[4], Ncoi, PM1[4], Nt1, PM2[4], Nt2;

  Stat = HDR[0];

  RunNum[1] = HDR[1];
  RunNum[0] = HDR[2];

  TRG[3] = HDR[3];
  TRG[2] = HDR[4];
  TRG[1] = HDR[5];
  TRG[0] = HDR[6];

  COI[3] = HDR[7];
  COI[2] = HDR[8];
  COI[1] = HDR[9];
  COI[0] = HDR[10];

  PM1[3] = HDR[11];
  PM1[2] = HDR[12];
  PM1[1] = HDR[13];
  PM1[0] = HDR[14];

  PM2[3] = HDR[15];
  PM2[2] = HDR[16];
  PM2[1] = HDR[17];
  PM2[0] = HDR[18];

  Run = RunNum[1]<<8 | RunNum[0];

  Nevt = TRG[3]<<24 | TRG[2]<<16 | TRG[1]<<8 | TRG[0];
  Ncoi = COI[3]<<24 | COI[2]<<16 | COI[1]<<8 | COI[0];
  Nt1  = PM1[3]<<24 | PM1[2]<<16 | PM1[1]<<8 | PM1[0];
  Nt2  = PM2[3]<<24 | PM2[2]<<16 | PM2[1]<<8 | PM2[0];

  printf("Run number %d Events %d Matches %d  PMT1 %d PMT2 %d\n", Run, Nevt, Ncoi, Nt1, Nt2);


/* copy the rest of the header into the buffer */

  result = fread (buffer, 1, (Hlength-46), fp);
  if (result != (Hlength-46)) {fputs ("Reading error", stderr); exit (3);}

  int GM2[(Hlength-46)];

  Value = buffer;

  for (i=0; i<(Hlength-46);i++) GM2[i] = (int)*Value++;

  for (i=0;i<(Hlength-46)/5;i++)
   {
     j= 5*i;

  /* number of the current Straw Active train */

       int Active = GM2[j];

  /* Clock cycle (after BOS) when Active train arrives */

       int Clock[4];

       Clock[3] = GM2[j+1];
       Clock[2] = GM2[j+2];
       Clock[1] = GM2[j+3];
       Clock[0] = GM2[j+4];

       int Time = Clock[3]<<24 | Clock[2]<<16 | Clock[1]<<8 | Clock[0];

//       printf("GM2 active %d  Clock count %d\n", Active, Time);
      }

/******************/
/* Loop over data */
/******************/

/* integer division to cut off corrupted tail */

  int Nblocks = (lSize - Hlength)/Blength;
  int Nloops = Nblocks / Nstations;

printf("Nblocks %d, blocks per station %d\n", Nblocks, Nloops);

  for (XYstat=0; XYstat<Nstations; XYstat++)
     {
       for (ind=0; ind<Nloops; ind++) 
          {
            int i;

	    int EDATA[Blength], MAC[6];

            int ADCX[155], ADCY[155], HX[16], HY[16];

/* copy the data portion into the buffer */

            result = fread (buffer, 1, Blength, fp);
            if (result != Blength) {fputs ("Reading error", stderr); exit (3);}

/* set pointer to the beginning of data field */

            Value = buffer;

            for (i=0; i<Blength; i++) EDATA[i] = (int)*Value++;

//if (XYstat==4 && ind==10) for(i=0;i<Blength;i++) printf("%02X\n", EDATA[i]);


/* MAC address associated with each Motherboard */

            for (i=0; i<6; i++) MAC[i] = EDATA[i];

/* XY station ID stored in the Master Card */

            int XYnum  = MAC[1]>>4 & 0xF;

/* Ethernet packet number to check readout integrity */

	    int PakNum = MAC[4]<<8 | MAC[5];


/* search through Data Block */

            i=6;
            int Goffset = Blength;

            while(i<(Blength-637))
              {
	        if (EDATA[i]==0xFF && EDATA[i+636]==0x55 && EDATA[i+637]==0xAA)
		  {
		    Goffset = i+1;
		    i = Blength;
		   }
		 i++;
		}

            if (Goffset==Blength) {}//printf("skip corrupt event %d\n", ind);
            else
               {

/* Amplitude Block */

                 int j=0;
                 for (i=0; i<620; i++)
                    {
                      int C(int n){return (EDATA[Goffset + n]>>7)&0x1;}

                     if (  (C(i+0)!=C(i+1)) &&
		           (C(i+1)==C(i+2)) &&
		           (C(i+2)==C(i+3)) &&
		           (C(i+3)==C(i+4)) )

		       { int A[4];

                        int I0 = (EDATA[i+Goffset+1]>>5)&0x3;
                        int I1 = (EDATA[i+Goffset+2]>>5)&0x3;
                        int I2 = (EDATA[i+Goffset+3]>>5)&0x3;
                        int I3 = (EDATA[i+Goffset+4]>>5)&0x3;

                         A[I0] = EDATA[i+Goffset+1]&0x1F;
                         A[I1] = EDATA[i+Goffset+2]&0x1F;
                         A[I2] = EDATA[i+Goffset+3]&0x1F;
                         A[I3] = EDATA[i+Goffset+4]&0x1F;

		         ADCX[j] = A[0]<<5 | A[1];
		         ADCY[j] = A[2]<<5 | A[3];
			 j++;

                       } /* end if C */

                    } /* end for Amplitudes */

// if ((XYstat==1)&&(ind==1)) for (i=0;i<j;i++) printf("%d %d\n", i, ADCX[i]);
// if ((XYstat==1)&&(ind==1)) for (i=0;i<j;i++) printf("%d %d\n", i, ADCY[i]);

/* Beetle header and strips */

                 double MX = ADCX[0]; /* Mean variable */
		 double MY = ADCY[0];

                 double RX = MX*MX;   /* RMS variable  */
		 double RY = MY*MY;

		 double SX = 0; /* StD variable */
		 double SY = 0;
		 
		 int iX = 0;    /* Pointer to X strips */
		 int iY = 0;    /* Pointer to Y strips */

                 int FX = 1;    /* Fuse to protect X pointer */
		 int FY = 1;    /* Fuse to protect Y pointer */

                 for (i=1; i<j; i++)
		    {
                      double n = (double)i+1;

		      MX = MX + ADCX[i];
		      MY = MY + ADCY[i];

                      RX = RX + ADCX[i]*ADCX[i];
		      RY = RY + ADCY[i]*ADCY[i];

                      SX = RX/n - MX*MX/n/n;
                      SY = RY/n - MY*MY/n/n;

/*

Acceptance:

X > Mean + N*sigma
X < Mean - N*sigma

(X-Mean)^2 > N^2*sigma^2  =>  (X-Mean)^2 > N^2 * (RMS^2 - Mean^2)

*/


                  /* hit the wall in X */

                      if ( (ADCX[i+1]-MX/n)*(ADCX[i+1]-MX/n) > Level*SX)
			   if (FX) {iX = i+1; FX = 0; /* Blow fuse X */} 

                  /* hit the wall in Y */

                      if ( (ADCY[i+1]-MY/n)*(ADCY[i+1]-MY/n) > Level*SY)
			   if (FY) {iY = i+1; FY = 0; /* Blow fuse Y */} 

		    } /* end for header and strips */

/* Header (pipeline counter etc) */

                 for (i=0;i<16;i++)
		      { HX[i] = ADCX[i+iX]; HY[i] = ADCY[i+iY]; }

/* Strips amplitudes */

                 for (i=0;i<(j-iX-16);i++) ADCX[i] = ADCX[i+iX+16];
                 for (i=0;i<(j-iY-16);i++) ADCY[i] = ADCY[i+iY+16];

/* TDC value (200ps resolution) */

                 int TDC = EDATA[Goffset+620]<<8 | EDATA[Goffset+621];

/* Temperature value */

                 int TMP = EDATA[Goffset+622]<<8 | EDATA[Goffset+623];

/* Run number */

                 int RUN = EDATA[Goffset+624]<<8 | EDATA[Goffset+625];

/* Event number */

                 int EVT = EDATA[Goffset+626]<<24 | (EDATA[Goffset+627]<<16) |
                           EDATA[Goffset+628]<<8  |  EDATA[Goffset+629];
/* Clock number */

                 int CLK = EDATA[Goffset+630]<<24 | (EDATA[Goffset+631]<<16) |
                           EDATA[Goffset+632]<<8  |  EDATA[Goffset+633];

/* Status register */

                 int STA = EDATA[Goffset+634];

/* XY station ID stored in each Motherboard */

                int XYadr = STA & 0xF;
		int ASIC0 = (STA>>4) & 0x1;
		int ASIC1 = (STA>>5) & 0x1;


/* Monitors */

// for (i=0; i<6; i++) printf(":%02X", EDATA[i]);
// printf("\n");

// if ((XYstat==1)&&(ind==1)) for (i=0;i<16;i++) printf("%d\n", HX[i]);
// if ((XYstat==1)&&(ind==1)) for (i=0;i<16;i++) printf("%d\n", HY[i]);

// if ((XYstat==1)&&(ind==1)) for (i=0;i<(j-iX-16);i++) printf("%d %d\n", i, ADCX[i]);
// if ((XYstat==1)&&(ind==1)) for (i=0;i<(j-iY-16);i++) printf("%d %d\n", i, ADCY[i]);

// printf("TDC  %d\n", TDC);
// printf("Temperature  %d\n", TMP);
// printf("Run   number %d\n", RUN);
// printf("Event number %d\n", EVT);
// printf("Clock number %d\n", CLK);
// printf("Status  %d\n", STA);

 printf("Station %d Status %02X, Event %d TDC %d\n", XYadr, STA, EVT, TDC);

               } /* end if Goffset */

          } /* end for Nevents */

    } /* end for XYstat */

  fclose (fp);
  free (buffer);
  return 0;
}