#include <stdlib.h>
#include <fcntl.h>
#include <stddef.h>
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <limits.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/param.h>
#include <sys/resource.h>
#include <sys/times.h>
#include <assert.h>
#include <strings.h>
#include <string.h>

#include "gdecomp.h"
#include "ctk.h"

#define MAXPOOL 20
#define DEBUG 0
#define FAKEID 0

FILE *TSlist;

long long int bytesRead = 0;

float gain[MAXDETNO], offset[MAXDETNO];
int maxevents;
long long int dtGT, dtEXT;
int minnGT, minnBGS, minnS800;

/*--------------------------------------------------------*/

void
CheckNoArgs (int required, int actual, char *str)
{

  if (required < actual)
    {
      printf ("argument problem with chat option\n");
      printf ("--> %s\n", str);
      printf ("required # arguments: %i\n", required - 1);
      printf ("actual   # arguments: %i\n", actual - 1);
      printf ("Please fix and try again, quitting...\n");
      exit (1);
    };

}

/*--------------------------------------------------------*/

int
readChatFile (char *name)
{

  /* declarations */

  FILE *fp, *fp1;
  char *p, *pc, str[STRLEN] = { '0' }, str1[STRLEN] =
  {
  '0'};
  char str2[STRLEN] = { '0' };
  int echo = 0, nret = 0, nni = 0, nn = 0;
  int i1, i;
  float r2, r3;
  int str_decomp (char *, int, int *);

  /* open chat file */

  if ((fp = fopen (name, "r")) == NULL)
    {
      printf ("error: could not open chat file: <%s>\n", name);
      exit (0);
    };
  printf ("chat file: <%s> open\n", name);
  printf ("\n");
  fflush (stdout);

  /* read content and act */

  pc = fgets (str, STRLEN, fp);

  while (pc != NULL)
    {
      if (echo)
        printf ("chat->%s", str);
      fflush (stdout);

      /* attemp to interpret the line */

      if ((p = strstr (str, "echo")) != NULL)
        {
          echo = 1;
          if (echo)
            printf ("will echo command lines\n");
          fflush (stdout);

        }
      else if (str[0] == 35)
        {
          /* '#' comment line, do nothing */
          nni--;                /* don't count as instruction */

        }
      else if (str[0] == 59)
        {
          /* ';' comment line, do nothing */
          nni--;                /* don't count as instruction */

        }
      else if (str[0] == 10)
        {
          /* empty line, do nothing */
          nni--;                /* don't count as instruction */

        }
      else if ((p = strstr (str, "maxevents")) != NULL)
        {
          nret = sscanf (str, "%s %i", str1, &maxevents);
          CheckNoArgs (nret, 2, str);
          printf ("will process a max of %i events\n", maxevents);
        }
      else if ((p = strstr (str, "dtGT")) != NULL)
        {
          nret = sscanf (str, "%s %lli", str1, &dtGT);
          CheckNoArgs (nret, 2, str);
          printf ("will use a GT time window of %lli\n", dtGT);
        }
      else if ((p = strstr (str, "dtEXT")) != NULL)
        {
          nret = sscanf (str, "%s %lli", str1, &dtEXT);
          CheckNoArgs (nret, 2, str);
          printf ("will use a BGS time window of %lli\n", dtEXT);
        }
      else if ((p = strstr (str, "minnGT")) != NULL)
        {
          nret = sscanf (str, "%s %i", str1, &minnGT);
          CheckNoArgs (nret, 2, str);
          printf ("minimum GT count before writeout is %i\n", minnGT);
        }
      else if ((p = strstr (str, "minnGS800")) != NULL)
        {
          nret = sscanf (str, "%s %i", str1, &minnS800);
          CheckNoArgs (nret, 2, str);
          printf ("minimum S800 count before writeout is %i\n", minnS800);
        }
      else if ((p = strstr (str, "minnBGS")) != NULL)
        {
          nret = sscanf (str, "%s %i", str1, &minnBGS);
          CheckNoArgs (nret, 2, str);
          printf ("minimum BGS count before writeout is %i\n", minnBGS);
        }
      else if ((p = strstr (str, "calfile")) != NULL)
        {
          nret = sscanf (str, "%s %s", str1, str2);
          CheckNoArgs (nret, 2, str);

          /* read calibration file */

          if ((fp1 = fopen (str2, "r")) == NULL)
            {
              printf ("could not open <%s>\n", str2);
            }
          else
            {
              printf ("%s is open\n", str2);

              fgets (str, 512, fp1);    /* header line */
              printf ("%s", str);

              i = fscanf (fp1, "\n%i %f %f", &i1, &r2, &r3);
              printf ("%i; %i %f %f\n", i, i1, r2, r3);
              nn = 0;
              while (i == 3 && i1 < MAXDETNO && i1 >= 0)
                {
                  offset[i1] = r2;
                  gain[i1] = r3;
                  nn++;
                  printf ("%i %f %f\n", i1, r2, r3);
                  i = fscanf (fp1, "\n%i %f %f", &i1, &r2, &r3);
                };
              fclose (fp1);
              printf ("read %i lines of ehi.cal\n", nn);

            };

        }
      else
        {

      /*-----------------------------*/
          /* chatscript read error point */
      /*-----------------------------*/

          printf ("line %2.2i in chat script, option :%s \n__not understood\n", nn, str);
          printf ("%i\n", str[0]);
          printf ("aborting\n");
          fflush (stdout);
          exit (0);
        };

      /* read next line in chat script */

      nn++;                     /* line counter */
      nni++;                    /* instruction counter */
      pc = fgets (str, STRLEN, fp);

    };



  /* done */

  fclose (fp);
  printf ("\n");
  printf ("chat file: <%s> closed\n", name);
  printf ("__processed %i sort instructions and %i lines\n", nni, nn);
  printf ("\n");
  fflush (stdout);
  return (0);

}

/*----------------------------------------------------------------------------*/

int
modData (CRYS_INTPTS * ev, GEBDATA * gd, float gain[MAXDETNO], float offset[MAXDETNO])
{

  /* declarations */

  int j, i1, i2, detno;

  if (gd->type != GEB_TYPE_DECOMP)
    {
//    printf("modData called, but data type was not GEB_TYPE_DECOMP, just return\n");
//    fflush(stdout);
      return (0);
    };

  /* apply gain/offset coefficients */

  i1 = (ev->crystal_id & 0x0003);
  i2 = ((ev->crystal_id & 0xfffc) >> 2);
  detno = i2 * 4 + i1;
  ev->tot_e *= gain[detno];
  ev->tot_e += offset[detno];

  for (j = 0; j < ev->num; j++)
    {
      ev->intpts[j].e *= gain[detno];
      ev->intpts[j].e += offset[detno];
    };

  /* done */

  return (0);

}


/*----------------------------------------------------------------------------*/

int
readHit (int in, CRYS_INTPTS * ev, GEBDATA * gd)
{

  /* declarations */

  int siz;
  int j, checksum, *iptr;
  static int nchecksum = 0;
  static int nn = 10, nr = 0, nts = 1000, nrr = 0;;
  static long long int lastTS = 0, lTS = 0;
  char str[127];

  assert (MAX_INTPTS == 16);

  if (nn > 0)
    printf ("\nreadHit: called, in=%p; nts=%i\n", (void *) in, nts);

  /* zero out */

  memset ((void *) ev, 0, sizeof (MAXDATASIZE));

  /* read event header in. It tells us the length */
  /* to use in the next read. For now we ignore the */
  /* type, but we may have to update that later */

  siz = read (in, (char *) gd, sizeof (GEBDATA));
  bytesRead += siz;
  if (siz != sizeof (GEBDATA))
    {
      printf ("readHit: cannot read GEBDATA, got siz=%i, not %i as expected\n", siz, sizeof (GEBDATA));
      return (2);
    };


  if (nn > 0)
    {
      printf ("readHit: read header of size %i, bytesRead=%lli, ", siz, bytesRead);
      printf ("gd->type/length/ts = %i, %i, %lli; ", gd->type, gd->length, gd->timestamp);
      printf ("dt= %lli\n", gd->timestamp - lastTS);
      if (gd->type == GEB_TYPE_DECOMP)
        lastTS = gd->timestamp;
      fflush (stdout);
    };


  if (nts > 0)
    {

      if ((int) (gd->timestamp - lTS) > 1000)
        fprintf (TSlist, "\n");
      if ((int) (gd->timestamp - lTS) < 0)
        fprintf (TSlist, "\n");

      get_GEB_Type_str (gd->type, str);
      fprintf (TSlist, "%s               , TS=%20lli, dt=%20lli\n", str, gd->timestamp, gd->timestamp - lTS);


      fflush (TSlist);

      if ((int) (gd->timestamp - lTS) > 1000)
        printf ("\n");
      if ((int) (gd->timestamp - lTS) < 0)
        printf ("\n");
      printf ("%s              , TS=%20lli, dt=%20lli\n", str, gd->timestamp, gd->timestamp - lTS);

      fflush (stdout);

      lTS = gd->timestamp;
    }
#if(0)
  else
    {
      printf ("debug stop\n");
      exit (0);
    };
#endif


  /* now read the data */

  if (nn > 0)
    {
      printf ("readHit: will attempt to read data body of length %i bytes(max=%i)\n", gd->length, MAXDATASIZE);
      fflush (stdout);
    };

  if (gd->length >= MAXDATASIZE)
    {
      printf ("gd->length= %i > %i (MAXDATASIZE)\n", gd->length, MAXDATASIZE);
      exit (1);
    };
//  assert (gd->length < MAXDATASIZE);
  siz = read (in, (char *) ev, gd->length);
  bytesRead += siz;
  if (nn > 0)
    printf ("readHit: read event of size %i\n", siz);
  if (siz != gd->length)
    {
      printf ("readHit: read error siz=%i != gd->length=%i\n", siz, gd->length);
      printf ("readHit: read %i events, bytesRead=%lli\n", nr, bytesRead);
      return (2);
    };

#if(0)

  /* check long GEB_TYPE_S800_RAW payloads */

  if (gd->type == GEB_TYPE_S800_RAW || 1)
//    if (gd->length>460)
    if (nchecksum < 10)
      {
        nchecksum++;
        checksum = 0;
        iptr = (int *) ev;
        for (j = 0; j < (gd->length / sizeof (int)); j++)
          {
            checksum += *iptr;
            printf ("%4i> %15i\n", j, *iptr);
            iptr++;
          };
        get_GEB_Type_str (gd->type, str);
        printf ("^^%s:: length= %6i, TS=%20lli, checksum=%10i\n", str, gd->length, gd->timestamp, checksum);
      };
#endif

  /* check, these are the only types we know about at this time */

//  assert ((gd->type == GEB_TYPE_DECOMP) || (gd->type == GEB_TYPE_BGS) || (gd->type == 8) );

  nr++;

#if(1)
  /* debug print what we read */

  if (nn > 0)
    {
      printf ("-----------------------------------------------------------------------\n");
      printf ("readHit: debug what we read, read no %i\n", nr);
      if (gd->type == GEB_TYPE_DECOMP)
        printCRYS_INTPTS (stdout, ev, gd);
      else if (gd->type == GEB_TYPE_BGS)
        printf ("BGS event with timestamp %lli\n", gd->timestamp);
      printf ("-----------------------------------------------------------------------\n");
    }
#endif

  /* done */

  if (nn > 0)
    nn--;
  if (nts > 0)
    nts--;

#if(0)
  if (nrr > 3)
    exit (0);
  nrr++;
#endif

  return (0);

}

/*----------------------------------------------------------------------------*/

int
main (int argc, char **argv)
{

  /* declarations */

  int i = 0, j = 0, nwritten = 0, nreads = 0, i1 = 0, i2, detno, st;
  long long int type_nbytes[100];
  int type_reads[100], nhit, nTSJumpBack = 0, nNotWritten = 0;
  int siz, nprint = 50, crystalhit[MAXDETNO];
  CRYS_INTPTS *ev[MAXPOOL];
  GEBDATA *gd[MAXPOOL];
  off_t inData = -1;
  off_t outData = -1;
  float dtbtevGT[16384], dtbtevEXT[16384];
  long long int deltaGT, deltaEXT;
  float r1;
  long long int TS_startEvent;
  int nBGS, nGT, nS800;
  int hitmat[20][20];
  double d1;
  char str[127];

  /* prototypes */

  int wr_spe (char *, int *, float *);

  if (argc != 4)
    {
      printf ("use %s <chatfile> <input> <output>\n", argv[0]);
      exit (0);
    };

  /*------------*/
  /* initialize */
  /*------------*/

  /* gain and offset coeff */

  for (i = 0; i < 20; i++)
    {
      type_reads[i] = 0;
      type_nbytes[i] = 0;
    };

  for (i = 0; i < MAXDETNO; i++)
    {
      gain[i] = 1;
      offset[i] = 0;
    };

  for (i = 0; i < 20; i++)
    for (j = 0; j < 20; j++)
      hitmat[i][j] = 0;

  /* read chat file with all support parameters */

  readChatFile (argv[1]);

  for (i = 0; i < MAXDETNO; i++)
    printf ("cal %3i: offset = %6.3f  gain = %9.4f\n", i, offset[i], gain[i]);

  /* help */

  if (argc == 1)
    {
      printf ("use: %s Merged.Mode2.dat 20000000000 40 400 Global_mod.dat ehi.cal \n", argv[0]);
      /*             0    1              2          3   4       5             6      */
      exit (0);
    };

  /* repeat pars */

  printf ("will construct %i events\n", maxevents);
  printf ("will use dtGT= %lli to define coincidence\n", dtGT);
  printf ("will use dtEXT= %lli to define coincidence\n", dtEXT);

  /* initialize */

  for (i = 0; i < MAXDETNO; i++)
    crystalhit[i] = 0;

  for (i = 0; i < 16384; i++)
    {
      dtbtevGT[i] = 0;
      dtbtevEXT[i] = 0;
    }

  for (i = 0; i < MAXPOOL; i++)
    {
      gd[i] = (GEBDATA *) calloc (1, sizeof (GEBDATA));
      ev[i] = (CRYS_INTPTS *) calloc (1, MAXDATASIZE);
    };

  /* open input file */

  inData = open (argv[2], O_RDONLY, 0);

  if (inData <= 0)
    {
      printf ("could not open input file %s\n", argv[2]);
      return (1);
    }
  else
    printf ("input file \"%s\" is open\n", argv[2]);

  /* open output file */

  outData = 0;
  /*                        + name of data file */
  /*                        |                   */
  outData = open ((char *) argv[3], O_WRONLY | O_CREAT | O_TRUNC, PMODE);
  if (outData == 0)
    {
      printf ("could not open output data file %s, quit!\n", argv[3]);
      exit (1);
    };

  printf ("output data file \"%s\" is open\n", argv[3]);
  fflush (stdout);

  /* simple ts list file */

  TSlist = fopen ("TS.list", "w");

  /* read first event */
  /* we requite this event to be GRETAINA data */
  /* assuming extenal data always comes later */

#if(0)
  /* very simple debug read */

  for (i = 0; i < 20; i++)
    st = readHit (inData, ev[0], gd[0]);

  if (1)
    exit (0);
#endif

  /* get the first event/prime the engine */

  st = readHit (inData, ev[0], gd[0]);
  if (st != 0)
    {
      printf ("could not read next event, quit\n");
      goto done;
    };
  nreads++;

  /* if decomposed gretina data, we can find detno etc... */


  if (gd[0]->type == GEB_TYPE_DECOMP)
    {
      modData (ev[0], gd[0], gain, offset);

      /* keep track of crystal hits */

      i1 = (ev[0]->crystal_id & 0x0003);
      i2 = ((ev[0]->crystal_id & 0xfffc) >> 2);
      detno = i2 * 4 + i1;
      if (detno < MAXDETNO)
        crystalhit[detno]++;

//          printCRYS_INTPTS (stdout,ev[0], gd[0]);
    }
  else
    {
      printf ("first data read was of type %i\n", gd[0]->type);
    };
  printf ("\n");

  TS_startEvent = gd[0]->timestamp;
  printf ("TS_startEvent = %lli\n", TS_startEvent);
  printf ("very first start data aquired\n\n");
  printf ("\n");

  /* now read events until timestamp is too far out */
  /* or we encounter an GEB_TYPE_S800PHYSDATA header */

  while (1)
    {

      j = 1;
      deltaGT = 0;
      deltaEXT = 0;
      while ((deltaGT < dtGT && deltaEXT < dtEXT) && j < MAXPOOL)
        {
          st = readHit (inData, ev[j], gd[j]);
          if (st != 0)
            {
              printf ("could not read event %i, finish up and quit\n", nreads);
              goto done;
            };


          /* check for timestamp that went back in time */

          /* Explanation: Especially for MSU data from Dirk */
          /* and Co, we get files that contain many runs. For */
          /* each run at MSU they reset the TS. Thus in the  */
          /* file there will be a lot of cases where the TS  */
          /* jumps back in time. We need to take care of that, */
          /* so we trap these cases here and manipulate  */
          /* TS_startEvent so the events are written out properly  */
          /* anyway. This is done in such a way that for files  */
          /* where the TS increases monotomicaly we still get  */
          /* the proper events out. */

          if (gd[j]->timestamp < TS_startEvent)
            {
              nTSJumpBack++;
              if (nTSJumpBack < 100)
                {
                  printf ("timestant went back in time!! ");
                  printf ("%i/%i\n", nwritten + 1, type_reads[9]);
                }
              else if (nTSJumpBack == 100)
                printf ("suspending TS jump warnings...\n");

              /* fix it, -1000 to make sure the current event is completed */

              TS_startEvent = gd[j]->timestamp - 1000;

//             exit(0);
            }

          /* statistics */

          nreads++;
          type_reads[gd[j]->type]++;
          type_nbytes[gd[j]->type] += gd[j]->length;

          /* if decomposed gretina data, we can find detno etc... */

          if (gd[j]->type == GEB_TYPE_DECOMP)
            {

              modData (ev[j], gd[j], gain, offset);

              /* keep track of crystal hits */

              i1 = (ev[j]->crystal_id & 0x0003);
              i2 = ((ev[j]->crystal_id & 0xfffc) >> 2);
              detno = i2 * 4 + i1;
              if (detno < MAXDETNO)
                crystalhit[detno]++;

              deltaGT = gd[j]->timestamp - TS_startEvent;
              if ((int) deltaGT < 16384 && (int) deltaGT >= 0)
                dtbtevGT[(int) deltaGT]++;

            }
          else if (gd[j]->type == GEB_TYPE_GT_MOD29 || gd[j]->type == GEB_TYPE_S800PHYSDATA)
            {
              deltaEXT = gd[j]->timestamp - TS_startEvent;
              if ((int) deltaEXT < 16384 && (int) deltaEXT >= 0)
                dtbtevEXT[(int) deltaEXT]++;
            }


#if(0)
          /* print event at this point */

          printf ("\n**current event (under construction) now looks like this:");
          printf ("j=%i\n", j);
          for (i = 0; i <= j; i++)
            {
              if (gd[i]->type == GEB_TYPE_DECOMP)
                printCRYS_INTPTS (stdout, ev[i], gd[i]);
              else if (gd[i]->type == GEB_TYPE_BGS)
                printf ("BGS event with timestamp %lli\n", gd[i]->timestamp);
            };
          printf ("deltaGT= %lli , deltaEXT= %lli\n", deltaGT, deltaEXT);
          printf ("TS_startEvent=%lli\n", TS_startEvent);
          printf ("**end of event\n");
#endif

          j++;


        };

      /* preserve start of next event */

//      TS_startEvent = gd[j - 1]->timestamp;
      TS_startEvent = gd[j - 1]->timestamp;
#if(0)
      printf ("setting start of next event to %lli\n", TS_startEvent);
#endif

      j--;



#if(0)
      printf ("... done\n");
#endif

      /* we now have an event */
      /* count how many GRETINA and EXTERNAL hits */
      /* we have in this event */

      nBGS = 0;
      nGT = 0;
      nS800 = 0;
      for (i = 0; i < j; i++)
        {
          if (gd[i]->type == GEB_TYPE_BGS)
            nBGS++;
          if (gd[i]->type == GEB_TYPE_DECOMP)
            nGT++;
          if (gd[i]->type == GEB_TYPE_S800PHYSDATA || gd[i]->type == GEB_TYPE_S800PHYSDATA)
            nS800++;
        };
#if(0)
      printf ("-------> nBGS=%i(%i),nGT=%i(%i),nS800=%i(%i)\n", nBGS, minnBGS, nGT, minnGT, nS800, minnS800);
#endif

      /* statistics */

      if (nGT < 20 && nBGS < 20)
        hitmat[nGT][nBGS]++;

      /* write out coincidence event if conditions are fulfilled */

      if (nBGS >= minnBGS && nGT >= minnGT && nS800 >= minnS800)
        {

          /* good event, write it out */

          if (nwritten < nprint)
            printf ("\n\n**writing coincidence event no %i with %i elements, type9=%i\n", nwritten + 1, j,
                    type_reads[9]);

          /* first we write the number of coincidence */
          /* events out that we have in this time bunch */

          siz = write (outData, (char *) &j, sizeof (int));

          /* next we write the data headers and actual data out */

          for (i = 0; i < j; i++)
            {

              siz = write (outData, (char *) gd[i], sizeof (GEBDATA));
              siz = write (outData, (char *) ev[i], gd[i]->length);

              /* debug printing */

              if (nwritten < nprint)
                {
                  if (gd[i]->type == GEB_TYPE_DECOMP)
                    printCRYS_INTPTS (stdout, ev[i], gd[i]);
                  else if (gd[i]->type == GEB_TYPE_DECOMP)
                    printf ("GEB_TYPE_DECOMP, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_RAW)
                    printf ("GEB_TYPE_RAW, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_TRACK)
                    printf ("GEB_TYPE_TRACK, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_BGS)
                    printf ("GEB_TYPE_BGS, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_S800_RAW)
                    printf ("GEB_TYPE_S800_RAW, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_NSCLnonevent)
                    printf ("GEB_TYPE_NSCLnonevent, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_GT_SCALER)
                    printf ("GEB_TYPE_GT_SCALER, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_GT_MOD29)
                    printf ("GEB_TYPE_GT_MOD29, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else if (gd[i]->type == GEB_TYPE_S800PHYSDATA)
                    printf ("GEB_TYPE_S800PHYSDATA, timestamp= %lli, length= %i\n", gd[i]->timestamp, gd[i]->length);
                  else
                    printf ("unknown data type= %i, timestamp= %lli, length= %i\n", gd[i]->type, gd[i]->timestamp,
                            gd[i]->length);
                };
            }
          if (nwritten < nprint)
            printf ("\n**done writing coincidence event no %i with %i elements\n\n", nwritten + 1, j);
          nwritten++;

//            if(1)exit(0);
        }
      else
        {
          if (nwritten < nprint)
            {
              printf ("THIS EVENT WAS NOT WRITTEN OUT\n");
              printf ("nBGS=%i ", nBGS);
              printf ("minnBGS=%i\n", minnBGS);
              printf ("nGT=%i ", nGT);
              printf ("minnGT=%i\n", minnGT);
              printf ("nS800=%i ", nS800);
              printf ("minnS800=%i\n", minnS800);

            };
          nNotWritten++;

        };

#if(0)
      /* report unusual conditions */

      if ((nS800 != 1) || (nGT == 0))
        {
          printf ("auch> nBGS=%i(%i),nGT=%i(%i),nS800=%i(%i)\n", nBGS, minnBGS, nGT, minnGT, nS800, minnS800);
//        exit(1);
        };
#endif

      /* BTW, are we done? */

      if (nwritten >= maxevents)
        goto done;

      /* reset, last event read is now first */

      memcpy (ev[0], ev[j], sizeof (CRYS_INTPTS));
      memcpy (gd[0], gd[j], sizeof (GEBDATA));

      TS_startEvent = gd[0]->timestamp;

    }

  /* done */

done:

  printf ("\n");
  printf ("crystal hitpattern:\n");
  nhit = 0;
  for (i = 0; i < MAXDETNO; i++)
    if (crystalhit[i] > 10)
      {
        printf ("%3i (module %2i, crystal %1i): %10i\n", i, i / 4, i - 4 * (i / 4), crystalhit[i]);
        nhit++;
      }
  printf (">>> seems %i crystals got hits\n", nhit);
  printf ("\n");

  printf ("\n\n ** wrote %7i events out from %7i reads (fraction %f)\n", nwritten, nreads, (float) nwritten / nreads);
  printf ("    %i events were not written out\n", nNotWritten);
  printf ("    fraction written coincidence events vs reads is %7.3f\n", (float) nreads / (float) nwritten);
  printf ("    wrote data to file: %s\n", argv[3]);
  printf ("we saw %i TSs jump back in time\n", nTSJumpBack);
  i1 = 16384;
  wr_spe ("dtbtevGT.spe", &i1, dtbtevGT);
  printf ("wrote dtbtevGT.spe\n");
  wr_spe ("dtbtevEXT.spe", &i1, dtbtevEXT);
  printf ("wrote dtbtevEXT.spe\n");
  printf ("wrote TS.list (for order debug)\n");

  d1 = (double) bytesRead / 1024000;
  printf ("read a total of %9.3f Mbytes from file\n", (float) d1);

#if(0)
  /* hitmat statistics */

  printf ("\n");
  printf ("GT vs BGS hit statistics\n");
  printf ("\n");

  i1 = 0;
  for (i = 0; i < 20; i++)
    for (j = 0; j < 20; j++)
      i1 += hitmat[i][j];

  i2 = 0;
  for (i = 0; i < 20; i++)
    for (j = 0; j < 20; j++)
      if (hitmat[i][j] > 10)
        {
          if (j > 0)
            i2 += i * hitmat[i][j];
          r1 = 100.0 * (float) hitmat[i][j] / i1;
          printf ("(nGT=%2i,nBGS%2i) %10i %9.3f %% %i\n", i, j, hitmat[i][j], r1, i2);
        };
  printf ("\n");
  printf ("we saw %10i gamma rays in coincidence\n", i2);
#endif

  printf ("any type events read  %10i (%9.2f %%)\n", nreads, 100.0);
  for (i = 0; i <= MAX_GEB_TYPE; i++)
    if (type_reads[i] > 0)
      {
        get_GEB_Type_str (i, str);
        printf ("geb type=%5i (%s): ", i, str);
        printf ("%10i reads ", type_reads[i]);
        r1 = 100.0 * (double) type_reads[i] / (double) nreads;
        printf ("(%9.2f %%) ", r1);
        printf ("and %9.3f MBytes", (float) ((double) type_nbytes[i] / 1024000));
        printf ("\n");
      };

  printf ("\ndone\n");
  return (0);


}