#include "stdlib.h"
#include "stdio.h"

int vert_checksum ( int data[273] );                 /* checksum calc */
int get_words ( FILE *fp, int nwords, int words[] ); /* read from input */
int ev_proc ( int data[] );                          /* process event   */
int find_word ( FILE *fp, int SearchFor );           /* read until SearchFor */
int hist[512][1024];

int idebug = 0;   /* increase this to reduce amount of output */
int NvParityError=0; /* counts vertical   parity errors */
int NhParityError=0; /* counts horizontal parity errors */

main (int argv, char *argc[])
{
/* This is a simple program written to try to decifer phenix raw  */
/* data format (prdf) files.    JPSullivan 18-Oct-99              */

/* modified format of event output to make it easier to understand*/
/* which output word is which.  JPSullivan  8-Dec-99              */

/* The calling sequence is    */
/* prdf filename.dat          */

/* *************************************************************** */

int nevent = 273; /* number of words in MCM output packet      */
int data[273];    /* to hold data from a single event          */

int Nev;          /* counts events read from file */

int istat;         /* status return from subroutine calls, <0 means error */

FILE *fp;            /* pointer to input prdf file name */
FILE *fhist;         /* pointer to histogram input prdf file name */

int i;               /* loop index */
int j;               /* loop index */

int FirstWord=0xFFFF;  /* first word in MVD packet is 16 bits all 1 */
/* **************************************************************** */

if (( fp = fopen( argc[1], "rb" )) == NULL )
  {
       if ( idebug<2 ) fprintf(stderr, "Error opening file, %s! \n", argc[1]);
  }
  else
  {
      /* file opened without errors, now start reading */
      if ( idebug<1 ) fprintf(stdout, "opened file %s\n", argc[1]);

      /* read events until we run out */

      for ( Nev=0; Nev<1000000; Nev++ )
      {
         /* get an event by looking for the expected first word */
         /* (low bit bits all 1) in an MVD data packet          */

         istat = find_word ( fp, FirstWord );
         if ( istat == -1 ) break;

         /* now the reset of the MCM data packet, 272 additional words  */

         data[0] = FirstWord;
         istat = get_words ( fp, nevent-1, &data[1] );
         /* if ( istat < 0 ) break; */

         istat = ev_proc ( data);    /* try to process the data */
         /* if ( istat < 0 ) break; */

      }
       
  }

  fhist = fopen("hist.hst", "w");

  for ( i=0 ; i<512; i++ )
  {
     for ( j=0; j<1023; j++ )
     {
        if ( hist[i][j]>0 ) {
           fprintf ( fhist, "%d %d %d\n", i,j,hist[i][j] );
        }
     }  
  }

  fprintf ( stdout,"\n histograms written to hist.hst\n");
  fprintf ( stdout,"\n total events read = %d\n",Nev);
  fprintf ( stdout,"Total vertical   parity errors seen = %d\n",NvParityError);
  fprintf ( stdout,"Total horizontal parity errors seen = %d\n",NhParityError);

  return;
}
int vert_checksum ( int data[273] )

/* This function was originally written by J.P.Sullivan 10-Jul-98.        */
/* This is a slight modification of the original. It is intended to       */
/* construct the checksum which appears in each MVD data packet. This     */
/* checksum is constructed from the data words and is summed across       */
/* all the words. The 1st bit of the checksum is the XOR of first bit     */
/* of all data words, the 2nd bit of the checksum from the 2nd bit of all */
/* data words, etc. This word is sometimes called the "vertical checksum" */
/* -- from which I take the name of the function.                         */

{
   int nbit;
   int nadc;
   int result;
   int exor;
   int ibits[16]={1,2,4,8,0x10,0x20,0x40,0x80,0x100,0x200,0x400,
                  0x800,0x1000,0x2000,0x4000,0x8000};
   int one=1;
  
   result=0;
   for (nbit=0; nbit<16; nbit++){
      exor=0;
      for ( nadc=0; nadc<271; nadc++){
         if ( ibits[nbit]&data[nadc] ) exor = exor^one;
      }
      if ( exor>0 ) result=result|ibits[nbit];
   }
   return result;
}
int get_words ( FILE *fp, int nwords, int words[] )
{
/* read nwords 32 bit words from the file fp and put the */
/* results in the array nwords.                          */
/* JPSullivan    18-Oct-1999                             */
/*********************************************************/

   int i;    /* loop over 32 bit words           */
   int ii;   /* i +- some offset where needed    */
   int j;    /* loop over bytes in a 32 bit word */

   int temp[4];   /* contains single bytes, used to assemble 32 bit words */

   int istat;     /* status return */
                  /* 0 means OK    */
                  /* -1 means EOF  */

   istat = 0;     /* default is OK */

   /* unpack header  */

   for ( i=0; i<6; i++) 
   {
      if      ( i==0 ) 
         { 
           if ( idebug<1 ) fprintf ( stdout, "\n\n Det ID  : ");
         }
      else if ( i==1 ) 
         {
           if ( idebug<1 ) fprintf ( stdout, " Evt num : ");
         }
      else if ( i==2 ) 
         {
           if ( idebug<1 ) fprintf ( stdout, " Mod add : ");
         }
      else if ( i==3 ) 
         {
           if ( idebug<1 ) fprintf ( stdout, "\n Flag wd : ");
         }
      else if ( i==4 ) 
         {
           if ( idebug<1 ) fprintf ( stdout, " N beam  : ");
         }
      else if ( i==5 ) 
         {
           if ( idebug<1 ) fprintf ( stdout, " AMUcells: ");
         }

      istat = get_one_word ( fp, i, &words[0] );

      if ( istat == -1 ) {
        if ( idebug < 2 ) fprintf ( stdout, " istat=-1, eof hit");
        return;
      }
   }

   /* unpack data  */

   for ( i=6; i<262; i++) 
   {
      ii = i-6;

      if ( (ii%8) == 0 ) { 
        if ( idebug<1) fprintf ( stdout, "\ndata: ");
      }

      istat = get_one_word ( fp, i, &words[0] );

      if ( istat == -1 ) {
        if ( idebug < 2 ) fprintf ( stdout, " istat=-1, eof hit");
        return;
      }
   }

   /* unpack trailing words   */

   for ( i=262; i<nwords; i++) 
   {
      ii = i-261;
      if ( ii==1 )     
         { 
            if ( idebug<1)  fprintf ( stdout, "\nuser: ");
         }
      else if ( ii==9 ) 
         {
            if ( idebug<1 ) fprintf ( stdout, "\nvpar: ");
         }
      else if ( ii==10) 
         {
            if ( idebug<1 ) fprintf ( stdout, "  zero:   ");
         }

      istat = get_one_word ( fp, i, &words[0] );

      if ( istat == -1 ) {
        if ( idebug < 2 ) fprintf ( stdout, " istat=-1, eof hit");
        return;
      }
   }

   return istat;
}
int get_one_word ( FILE *fp, int index, int words[] )
{
/* read one 32 bit word from the file fp and put the     */
/* results in the array words at position index.         */
/* JPSullivan    8-Dec-1999                              */
/*********************************************************/

   int j;    /* loop over bytes in a 32 bit word */

   int temp[4];   /* contains single bytes, used to assemble 32 bit words */

   int istat;     /* status return */
                  /* 0 means OK    */
                  /* -1 means EOF  */

   istat = 0;     /* default is OK */

   for ( j=0; j<4; j++)
   {
     temp[j] = getc(fp);
     if ( temp[j] == EOF ) {
       istat=-1; 
       return;
     }
   }
   words[index] = (temp[0] << 24) | (temp[1] << 16) |
                  (temp[2] <<  8) | temp[3];
   if ( idebug < 1 ) fprintf ( stdout, "%9X ",words[index]);

   return istat;
}
int ev_proc ( int data[] )
{

/* Simple routine to process MCM output packets */
/* written by JPSullivan 18-Oct-1999            */
/************************************************/

   int v_calc;        /* vertical parity calculated */
   int v_parity;      /* vertical parity from data  */
   int mask16bits = 0xFFFF;  /* to mask off lowest 16 bits */

   int istat;       /*  status return */
                    /*  0 means OK    */
                    /* +1 means v parity error */
                    /* +2 means h parity error */
                    /* -1 first word wrong   */
                    /* -2 wrong detector ID  */

   int all_one_word=0;         /* element number in array */
   int all_one_value = 0xFFFF; /* expected value          */
   int det_id_word=1;          /* element number in array */
   int det_id_value=2;         /* expected value          */

   int NuserWord1=263;         /* index to 1st  user word in data */
   int NuserWord8=270;         /* index to last user word in data */

   int maskbit10 = 0x400;      /* bit 10 set */
   int maskbit11 = 0x800;      /* bit 11 set */
   int maskbit12 = 0x1000;     /* bit 12 set */
   int maskbit13 = 0x2000;     /* bit 13 set */
   int maskbit14 = 0x4000;     /* bit 14 set */

   int i;                      /* loop index */

   int NhistOffset;   /* 0 for even MCMs(0,2,4), 256 for odd (1,3,5) */
   int ADCvalue;      /* single ADC value extracted from data array  */
   int NADC1=7;       /* pointer to first ADC value in data array    */

/***************************************************/

   istat = 0;

   /* check that a few words are in the right place */

   if ( (data[all_one_word]&mask16bits) != all_one_value ) {
      if ( idebug < 2 ) fprintf ( stdout, 
                "\nfirst word should be %x, is actually=%x",
                all_one_value, data[all_one_word]&mask16bits);
      istat = -1;
      return istat;
   }

   if ( (data[det_id_word]&mask16bits) != det_id_value ) {
      if ( idebug < 2 ) fprintf ( stdout, 
                        "\nDetector ID should be %x, is actually=%x",
                        det_id_value, data[det_id_word]&mask16bits);
      istat = -2;
      return istat;
   }   

   if ( (data[NuserWord8]&maskbit10) == 0 ) {
      if ( idebug < 2 ) fprintf ( stdout,
                       "\nValid Data bit is not set");
      istat = -3;
      return istat;
   }

   if ( (data[NuserWord8]&maskbit11) == 1 ) {
      if ( idebug < 2 ) fprintf ( stdout,
                       "\npacket too short");
      istat = -4;
      return istat;
   }

   if ( (data[NuserWord8]&maskbit12) == 1 ) {
      if ( idebug < 2 ) fprintf ( stdout,
                       "\nno stop sequence");
      istat = -5;
      return istat;
   }

   if ( (data[NuserWord8]&maskbit13) == 1 ) {
      if ( idebug < 2 ) fprintf ( stdout,
                       "\ntwo consecutive start sequences");
      istat = -6;
      return istat;
   }

   /* the packet format looks OK, so check parity */

   /* vertical parity */

   v_parity = data[271] & mask16bits;
   v_calc   = vert_checksum ( data );

   if ( v_parity == v_calc ) {
      if ( idebug < 1 ) fprintf ( stdout, "\n vertical parity word OK");
   }
      else
   {
      if ( idebug < 2 ) fprintf ( stdout, 
                       "\n vertical parity ERROR, data=%X, calc=%X",
                       v_parity, v_calc);
      istat+=1;
      NvParityError++;
   }

   /* horizontal parity */

   for ( i=0; i<7 ; i++ )
   {
      if ( (data[NuserWord1+i] & mask16bits) != 0 ) 
      {
        NhParityError++;
        istat +=2;
        if ( idebug < 2 ) fprintf ( stdout, 
                          "\n horizontal parity ERROR, word=%x",
                          data[NuserWord1+i]);
        break;
      }
   }

   if ( istat==0 ) 
   {
      /* if everything seems OK, put the data into histograms */
     
      if ( (data[NuserWord1+i] & maskbit14) ) NhistOffset = 256;
      else                                    NhistOffset = 0;

      for ( i=0; i<256; i++ )
      {
         ADCvalue = data[NADC1+i] & mask16bits;
         if ( (ADCvalue >= 0) && (ADCvalue <=1023) )
         {
            hist[NhistOffset+i][ADCvalue]++;     /* increment histogram */
         }
      }


      
   }

   return istat;
}
int find_word ( FILE *fp, int SearchFor )
{
/* read one 32 bit word at a time from file fp, until   */
/* SearchFor is found.                                  */
/* JPSullivan    22-Oct-1999                            */
/********************************************************/

   int i;    /* loop over 32 bit words           */
   int j;    /* loop over bytes in a 32 bit word */

   int temp[4];   /* contains single bytes, used to assemble 32 bit words */
   int DataValue; /* contains 32 bit word assembled from 4 bytes */

   int istat;     /* status return */
                  /* 0 means OK    */
                  /* -1 means EOF  */

   istat = 0;     /* default is OK */

   while ( (SearchFor != DataValue) && (istat != -1) ) 
   {
      for ( j=0; j<4; j++)
      {
        temp[j] = getc(fp);
        if ( temp[j] == EOF ) {
          istat=-1; 
          if ( idebug < 2 ) fprintf ( stdout, " find_word: istat=-1, EOF hit");
          break;
        }
      }
      DataValue = (temp[0] << 24) | (temp[1] << 16) |
                  (temp[2] <<  8) | temp[3];
   }

   return istat;
}