#include <stdlib.h>
#include <fcntl.h>
#include <stddef.h>
#include <stdio.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.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 <string.h>

#include "ctk.h"

#define DEBUG 0
#define ZEROSUPPRESS 1

/* globals */

ANAPAR Pars;
//CLUSTER_INTPTS Clstr[MAXCLUSTERHITS];
//int nClusters;
float rad2deg;
int dethit[MAXDETNO];

off_t infile;
double mean_pol, mean_azi;
int mean_n;


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

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
unpackBGS (TRACK_STRUCT * track, BGS_STRUCT * bgs)
{

  /* declarations */

  int k, i, j;
  PAYLOAD *ptinp;
  GEBDATA *ptgd;
  unsigned short int *i2;
  struct bgsp
  {
    unsigned short int par;
    unsigned short int val;
  };

  /* init  */

  bgs->nEsi = 0;

  /* find BGS data */

  ptinp = track->payload;
  ptgd = track->gd;
  for (k = 0; k < track->n; k++)
    {
      if (ptgd->type == GEB_TYPE_BGS)
        {
//        printf("found BGS data of length %i, bytes or %i short ints\n", ptgd->length, ptgd->length/2);

          i2 = (unsigned short int *) ptinp;
//        printf("    --- %i\n", *i2);

#if(0)
          /* just print and quit (DEBUG) */

          j = ptgd->length / sizeof (unsigned short int);
          for (i = 0; i < j; i++)
            {
              printf ("%3i: %8i\n", i, *i2);
              i2++;
            };
          if (1)
            exit (0);
#endif

          /* extract Si energies */

          j = ptgd->length / (2 * sizeof (unsigned short int));

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

//           printf("par=%8i-> val=%8i\n",*i2, *(i2+1));

              /* extract BGS energy per I-Y code */

              bgs->Esi[bgs->nEsi] = -1;

              if (*i2 >= 64 && *i2 <= 95)
                bgs->Esi[bgs->nEsi] = *(i2 + 1);
              else if (*i2 >= 128 && *i2 <= 159)
                bgs->Esi[bgs->nEsi] = *(i2 + 1);
              else if (*i2 >= 192 && *i2 <= 223)
                bgs->Esi[bgs->nEsi] = *(i2 + 1);
              else if (*i2 >= 256 && *i2 <= 287)
                bgs->Esi[bgs->nEsi] = *(i2 + 1);

              if (bgs->Esi[bgs->nEsi] != -1)
                {
//                  printf("bgs->Esi[%i]=%8i (par=%i)\n",bgs->nEsi,bgs->Esi[bgs->nEsi], *i2);
                  bgs->nEsi++;
                };

              /* next */

              i2 += 2;

            };
        };

      /* next */

      ptinp++;
      ptgd++;

    };

  /* done */

  return (0);

}

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


int
ProcessExternal (TRACK_STRUCT * track, int EXT_tac, BGS_STRUCT * bgs, int *ntac, double tac[100])
{

  /* in this routine we can fish out ProcessExternal data hidden in the track */
  /* structure and make a determination of whether to call ctkROOTBin */
  /* or not. Return '1' to bin, '0' to not bin */

  /* declarations */

  int ok = 0, k;

  /* quit here if we are not asked to process external data */

  if (Pars.gateexternal != 1)
    return (1);

  /*------------------------*/
  /* GT time tac requirement */
  /*------------------------*/
  /* not done here anymore */



//  printf("EXT_tac=%i\n", EXT_tac);
//  if (EXT_tac >= Pars.EXT_tac_lo && EXT_tac <= Pars.EXT_tac_hi)
//    ok = 1;

//  if (!ok)
//    return (0);

  ok = 1;


  /*--------------------------------------------*/
  /* before we go further: unpack the BGS data */
  /*--------------------------------------------*/

//  unpackBGS (track, bgs);


  /*----------------------*/
  /* BGS Esi requirements */
  /*----------------------*/

  ok = 0;
  for (k = 0; k < bgs->nEsi; k++)
    {
      if (bgs->Esi[k] > Pars.EXT_e_lo && bgs->Esi[k] <= Pars.EXT_e_hi)
        ok++;
    };
//  printf("ok=%i,bgs->nEsi=%i \n", ok,bgs->nEsi);

  /* for now require all or just one energies to be in range */
  /* what is that the correct thing to do ??? */
  /* :: require just one in the proper range */

//  if (ok == bgs->nEsi)
  if (ok > 0)
    ok = 1;
  else
    ok = 0;

  /* return evaluation */

  if (ok)
    return (1);
  else
    return (0);

};

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

int
findDetectorNumbers (TRACK_STRUCT * track)
{

  /* declarations */

  PAYLOAD *ptinp;
  GEBDATA *ptgd;
  CRYS_INTPTS *gtinp;
  int i, crystalno, moduleno, i1;

  /* mod data */

  ptinp = track->payload;
  ptgd = track->gd;
  for (i = 0; i < track->n; i++)
    {
      gtinp = (CRYS_INTPTS *) ptinp;


      if (ptgd->type == GEB_TYPE_DECOMP)
        {
          gtinp = (CRYS_INTPTS *) ptinp;

          crystalno = (gtinp->crystal_id & 0x0003);
          moduleno = (gtinp->crystal_id & 0xfffc) >> 2;
          i1 = moduleno * 4 + crystalno + 1;
          if (i1 >= 0 && i1 < MAXDETNO)
            {
              dethit[i1]++;
              Pars.detNo[i] = i1;
            };
        };
      ptgd++;
      ptinp++;

    };

  /* done */

  return (0);

}


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

int
readChatFile (char *name)
{

  /* declarations */

  FILE *fp;
  char *p, *pc, str[STRLEN] = { '0' }, str1[STRLEN] =
  {
  '0'};
  char str2[STRLEN] = { '0' };
  int echo = 0, nret = 0, nni = 0, nn = 0;
  int i1, j;
  float r1, r2, rr;
  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, "infile")) != NULL)
        {
          nret = sscanf (str, "%s %s", str1, str2);
          CheckNoArgs (nret, 2, str);

          /* open the output file */

          infile = open (str2, O_RDONLY, 0);

          if (infile <= 0)
            {
              printf ("could not open output file %s\n", str2);
              return (1);
            }
          else
            printf ("input file \"%s\" is open\n", str2);

        }
      else if ((p = strstr (str, "rmin")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.rmin);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "rmax")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.rmax);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "rootfile")) != NULL)
        {
          nret = sscanf (str, "%s %s", str1, Pars.rootFileName);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "s800_pid_x")) != NULL)
        {
          nret = sscanf (str, "%s %i %f %f", str1, &Pars.S800_PID_NX, &Pars.S800_PID_XLO, &Pars.S800_PID_XHI);
          CheckNoArgs (nret, 4, str);
        }
      else if ((p = strstr (str, "s800_pid_y")) != NULL)
        {
          nret = sscanf (str, "%s %i %f %f", str1, &Pars.S800_PID_NY, &Pars.S800_PID_YLO, &Pars.S800_PID_YHI);
          CheckNoArgs (nret, 4, str);
        }
      else if ((p = strstr (str, "ytascale")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.ytascale);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "greta_ata")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.greta_ata);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "greta_bta")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.greta_bta);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "target_x")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.target_x);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "target_y")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.target_y);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "target_z")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.target_z);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "minazi")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.minazi);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "maxazi")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.maxazi);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "minpol")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.minpol);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "maxpol")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.maxpol);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "PIDwin")) != NULL)
        {
          nret = sscanf (str, "%s %s", str1, Pars.PIDwinfile);
          CheckNoArgs (nret, 2, str);
          Pars.havePIDwin = 1;
        }
      else if ((p = strstr (str, "simpleDopplerCor")) != NULL)
        {
          CheckNoArgs (nret, 1, str);
          Pars.simpleDopplerCor = 1;
        }
      else if ((p = strstr (str, "exit")) != NULL)
        {
          fclose (fp);
          return (0);
        }
      else if ((p = strstr (str, "enabled")) != NULL)
        {
          nret = sscanf (str, "%s %s", str1, str2);
          CheckNoArgs (nret, 2, str);
          str_decomp (str2, MAXDETNO + 1, Pars.enabled);
          for (j = 0; j < MAXDETNO; j++)
            if (!Pars.enabled[j])
              printf ("detector %3i is disabled\n", j);
            else
              printf ("detector %3i is enabled\n", j);
        }
      else if ((p = strstr (str, "nprint")) != NULL)
        {
          nret = sscanf (str, "%s %i", str1, &Pars.nprint);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "requireext")) != NULL)
        {
          Pars.gateexternal = 1;
          printf ("will gate with external data, %i\n", Pars.gateexternal);
        }
      else if ((p = strstr (str, "exttaclim")) != NULL)
        {
          nret = sscanf (str, "%s %i %i", str1, &Pars.EXT_tac_lo, &Pars.EXT_tac_hi);
          CheckNoArgs (nret, 3, str);
          printf ("will sort with BGS tac time limits: from %i to %i\n", Pars.EXT_tac_lo, Pars.EXT_tac_hi);
        }
      else if ((p = strstr (str, "bgsesilim")) != NULL)
        {
          nret = sscanf (str, "%s %i %i", str1, &Pars.EXT_e_lo, &Pars.EXT_e_hi);
          CheckNoArgs (nret, 3, str);
          printf ("will sort with BGS silicon E limits: from %i to %i\n", Pars.EXT_e_lo, Pars.EXT_e_hi);
        }
      else if ((p = strstr (str, "beamdir")) != NULL)
        {
          nret = sscanf (str, "%s %f %f %f", str1, &Pars.beamdir[0], &Pars.beamdir[1], &Pars.beamdir[2]);
          CheckNoArgs (nret, 4, str);
          rr =
            Pars.beamdir[0] * Pars.beamdir[0] + Pars.beamdir[1] * Pars.beamdir[1] + Pars.beamdir[2] * Pars.beamdir[2];
          rr = sqrtf (rr);
          printf ("will use beam direction %f %f %f\n", Pars.beamdir[0], Pars.beamdir[1], Pars.beamdir[2]);
          printf ("vector length= %f, normalizing\n", rr);
          Pars.beamdir[0] /= rr;
          Pars.beamdir[1] /= rr;
          Pars.beamdir[2] /= rr;
          printf ("will use beam direction %f %f %f\n", Pars.beamdir[0], Pars.beamdir[1], Pars.beamdir[2]);
        }
      else if ((p = strstr (str, "beta")) != NULL)
        {
          nret = sscanf (str, "%s %f", str1, &Pars.beta);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "fomlimits")) != NULL)
        {
          nret = sscanf (str, "%s %i %f %f", str1, &i1, &r1, &r2);
          CheckNoArgs (nret, 3, str);
          Pars.fomlo[i1] = r1;
          Pars.fomhi[i1] = r2;
        }
      else if ((p = strstr (str, "ndetlimits")) != NULL)
        {
          nret = sscanf (str, "%s %f %f", str1, &Pars.ndetlimlo, &Pars.ndetlimhi);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "binch2d")) != NULL)
        {
          nret = sscanf (str, "%s %i", str1, &Pars.nBinCh2d);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "SMAPsiz")) != NULL)
        {
          nret = sscanf (str, "%s %i", str1, &Pars.nSMAP);
          CheckNoArgs (nret, 2, str);
        }
      else if ((p = strstr (str, "testdata")) != NULL)
        {
          nret = sscanf (str, "%s %s", str1, Pars.testdatafn);
          CheckNoArgs (nret, 2, str);
          Pars.havetestdata = 1;
          printf ("will load \"%s\" as test data\n", Pars.testdatafn);
        }
      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
main (int argc, char **argv)
{
  /* declarations */

  int i, j, k, HaveChatFile = 0, st = 0;
  int EXT_tac;
  float polang[MAXCLUSTERHITS], rr;
  float doppler_factor[MAXCLUSTERHITS];
  long long int meanGT_time, meanBGS_time, meanS800_time, nmeanGT_time;
  int nmeanBGS_time, nmeanS800_time;
  float r1;
  char *p, str[132];
  char ChatFileName[STRLEN];
  TRACK_STRUCT track;
  BGS_STRUCT bgs;
  float pi;
  FILE *fp;
  PAYLOAD *ptinp;
  GEBDATA *ptgd;
  int getTrackedEvent (TRACK_STRUCT *);
  int hitmat[20][20], i1, i2;

  /* inform of version of root used */

  printf ("from %s on %s\n", __FILE__, __DATE__);

  /* initialize */

  Pars.havetestdata = 0;
  pi = 4.0 * atan (1.0);
  rad2deg = pi / 180;

  Pars.GGMAX = 2000;
  mean_pol = 0;
  mean_azi = 0;
  mean_n = 0;
  Pars.rmin = 0;
  Pars.rmax = 30;
  Pars.minazi = 0;
  Pars.maxazi = 360;
  Pars.minpol = 0;
  Pars.maxpol = 180;


  /* default PID limits */

  Pars.S800_PID_NX = 512;
  Pars.S800_PID_XLO = 0;
  Pars.S800_PID_XHI = 3000;

  Pars.S800_PID_NY = 512;
  Pars.S800_PID_YLO = 0;
  Pars.S800_PID_YHI = 3000;

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

  Pars.ytascale = 1;

#if(0)

  /* THIS if for ER setup */

  /* crystal polar angles per documentation */
  /* eventually need all of them in here.... */

  for (i = 1; i <= MAXDETNO; i++)
    Pars.detpolang[i] = 90;

  /* module # 9 */

  Pars.detpolang[36] = 67.71;
  Pars.detpolang[37] = 55.02;
  Pars.detpolang[38] = 48.92;
  Pars.detpolang[39] = 64.18;

  /* module # 17 */

  Pars.detpolang[68] = 87.15;
  Pars.detpolang[69] = 105.79;
  Pars.detpolang[70] = 92.85;
  Pars.detpolang[71] = 74.21;

  /* module # 19 */

  Pars.detpolang[76] = 87.15;
  Pars.detpolang[77] = 105.79;
  Pars.detpolang[78] = 92.85;
  Pars.detpolang[79] = 74.21;

  /* module # 20 */

  Pars.detpolang[80] = 97.01;
  Pars.detpolang[81] = 100.78;
  Pars.detpolang[82] = 82.99;
  Pars.detpolang[83] = 79.22;

  /* module # 24 */

  Pars.detpolang[96] = 112.29;
  Pars.detpolang[97] = 124.98;
  Pars.detpolang[98] = 131.08;
  Pars.detpolang[99] = 115.82;

  /* module # 25 */

  Pars.detpolang[100] = 112.99;
  Pars.detpolang[101] = 124.98;
  Pars.detpolang[102] = 131.08;
  Pars.detpolang[103] = 115.82;

  /* module # 28 */

  Pars.detpolang[112] = 144.7;
  Pars.detpolang[113] = 132.4;
  Pars.detpolang[114] = 149.4;
  Pars.detpolang[115] = 163.35;

  /* module # 29 */

  Pars.detpolang[116] = 144.7;
  Pars.detpolang[117] = 132.4;
  Pars.detpolang[118] = 149.4;
  Pars.detpolang[119] = 163.35;

#endif

#if(1)

  /* this is for BGS setup */

  Pars.detpolang[4] = 92.853203;
  Pars.detpolang[5] = 74.206596;
  Pars.detpolang[6] = 87.146698;
  Pars.detpolang[7] = 105.792999;
  Pars.detpolang[8] = 82.993500;
  Pars.detpolang[9] = 79.221397;
  Pars.detpolang[10] = 97.005798;
  Pars.detpolang[11] = 100.778000;
  Pars.detpolang[12] = 112.286003;
  Pars.detpolang[13] = 124.982002;
  Pars.detpolang[14] = 131.082001;
  Pars.detpolang[15] = 115.819000;
  Pars.detpolang[16] = 149.401993;
  Pars.detpolang[17] = 163.350006;
  Pars.detpolang[18] = 144.699997;
  Pars.detpolang[19] = 132.397995;
  Pars.detpolang[20] = 131.082001;
  Pars.detpolang[21] = 115.819000;
  Pars.detpolang[22] = 112.287003;
  Pars.detpolang[23] = 124.982002;
  Pars.detpolang[24] = 48.917599;
  Pars.detpolang[25] = 64.180603;
  Pars.detpolang[26] = 67.713203;
  Pars.detpolang[27] = 55.017601;
  Pars.detpolang[28] = 87.146301;
  Pars.detpolang[29] = 105.792999;
  Pars.detpolang[30] = 92.852600;
  Pars.detpolang[31] = 74.206100;
  Pars.detpolang[32] = 149.401993;
  Pars.detpolang[33] = 163.348999;
  Pars.detpolang[34] = 144.699005;
  Pars.detpolang[35] = 132.397995;
  Pars.detpolang[36] = 144.699997;
  Pars.detpolang[37] = 132.399002;
  Pars.detpolang[38] = 149.401993;
  Pars.detpolang[39] = 163.350006;
  Pars.detpolang[40] = 82.994102;
  Pars.detpolang[41] = 79.222000;
  Pars.detpolang[42] = 97.006401;
  Pars.detpolang[43] = 100.778999;
  Pars.detpolang[44] = 67.713799;
  Pars.detpolang[45] = 55.018101;
  Pars.detpolang[46] = 48.918201;
  Pars.detpolang[47] = 64.181297;
  Pars.detpolang[48] = 30.598200;
  Pars.detpolang[49] = 16.649900;
  Pars.detpolang[50] = 35.299900;
  Pars.detpolang[51] = 47.601398;
  Pars.detpolang[52] = 30.597601;
  Pars.detpolang[53] = 16.649500;
  Pars.detpolang[54] = 35.299500;
  Pars.detpolang[55] = 47.600800;
  Pars.detpolang[56] = 48.917301;
  Pars.detpolang[57] = 64.180496;
  Pars.detpolang[58] = 67.712898;
  Pars.detpolang[59] = 55.017101;
  Pars.detpolang[60] = 82.993301;
  Pars.detpolang[61] = 79.221397;
  Pars.detpolang[62] = 97.005699;
  Pars.detpolang[63] = 100.778000;
  Pars.detpolang[64] = 112.286003;
  Pars.detpolang[65] = 124.982002;
  Pars.detpolang[66] = 131.082001;
  Pars.detpolang[67] = 115.819000;
  Pars.detpolang[68] = 144.699997;
  Pars.detpolang[69] = 132.399002;
  Pars.detpolang[70] = 149.401993;
  Pars.detpolang[71] = 163.350006;
  Pars.detpolang[72] = 144.699997;
  Pars.detpolang[73] = 132.399002;
  Pars.detpolang[74] = 149.403000;
  Pars.detpolang[75] = 163.350006;
  Pars.detpolang[76] = 131.082993;
  Pars.detpolang[77] = 115.820000;
  Pars.detpolang[78] = 112.287003;
  Pars.detpolang[79] = 124.983002;
  Pars.detpolang[80] = 92.853600;
  Pars.detpolang[81] = 74.207100;
  Pars.detpolang[82] = 87.147301;
  Pars.detpolang[83] = 105.793999;
  Pars.detpolang[84] = 35.300499;
  Pars.detpolang[85] = 47.601799;
  Pars.detpolang[86] = 30.598499;
  Pars.detpolang[87] = 16.650499;
  Pars.detpolang[88] = 30.597500;
  Pars.detpolang[89] = 16.649700;
  Pars.detpolang[90] = 35.299702;
  Pars.detpolang[91] = 47.600800;
  Pars.detpolang[92] = 92.852798;
  Pars.detpolang[93] = 74.206299;
  Pars.detpolang[94] = 87.146698;
  Pars.detpolang[95] = 105.792999;
  Pars.detpolang[96] = 131.082001;
  Pars.detpolang[97] = 115.819000;
  Pars.detpolang[98] = 112.287003;
  Pars.detpolang[99] = 124.983002;
  Pars.detpolang[100] = 97.006699;
  Pars.detpolang[101] = 100.778999;
  Pars.detpolang[102] = 82.994301;
  Pars.detpolang[103] = 79.222397;
  Pars.detpolang[104] = 30.598101;
  Pars.detpolang[105] = 16.650400;
  Pars.detpolang[106] = 35.300400;
  Pars.detpolang[107] = 47.601398;
  Pars.detpolang[108] = 67.713303;
  Pars.detpolang[109] = 55.017399;
  Pars.detpolang[110] = 48.917801;
  Pars.detpolang[111] = 64.181099;
  Pars.detpolang[112] = 97.006302;
  Pars.detpolang[113] = 100.778000;
  Pars.detpolang[114] = 82.993797;
  Pars.detpolang[115] = 79.222000;
  Pars.detpolang[116] = 87.147301;
  Pars.detpolang[117] = 105.793999;
  Pars.detpolang[118] = 92.853500;
  Pars.detpolang[119] = 74.207001;
  Pars.detpolang[120] = 48.918301;
  Pars.detpolang[121] = 64.181503;
  Pars.detpolang[122] = 67.713898;
  Pars.detpolang[123] = 55.018002;

#endif

  /* help ? */

  if (argc < 2)
    {
      printf ("use: %s -chat file [-version] [-help]\n", argv[0]);
      exit (0);
    };

  /* now read the chat script file */

  j = 1;
  if (argc > 1)
    while (j < argc)
      {

        if ((p = strstr (argv[j], "-version")) != NULL)
          {
            printf ("not sure, now using svn...\n");
            exit (0);
          }
        else if ((p = strstr (argv[j], "-help")) != NULL)
          {
            printf ("\n");
            printf ("ctk will be documented on the WWW, URL:\n");
            printf ("\n");
            printf ("  http://www.phy.anl.gov/somewhere_to_be_done\n");
            printf ("\n");
            exit (0);
          }
        else if ((p = strstr (argv[j], "-chat")) != NULL)
          {
            j++;
            strcpy (ChatFileName, argv[j++]);
            printf ("will read analysis instructions from chatfile: %s\n", ChatFileName);
            HaveChatFile = 1;
          }
        else
          {
            printf ("command line argument not understood!\n");

            printf ("%s: I was called as: \n--->[", argv[0]);
            for (i = 0; i < argc; i++)
              {
                printf ("%s ", argv[i]);
                fflush (stdout);
              }
            exit (0);

          }
      };

  /* read chat file - or die */

  if (HaveChatFile)
    readChatFile (ChatFileName);
  else
    {
      printf ("error: you must specify a chat script\n");
      exit (1);
    };


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

  for (i = 1; i <= MAXDETNO; i++)
    {
      /* these are the simple per detector */
      /* doppler correction factors; */
      /* not what we use when we have done tracking */

      rr = 1.0 - Pars.beta * Pars.beta;
      Pars.detdopfac[i] = sqrt ((double) rr) / (1.0 - Pars.beta * cos ((double) Pars.detpolang[i] * rad2deg));
    };

  for (i = 1; i <= MAXDETNO; i++)
    {
      printf ("det %3i, pol= %9.2f, doppfactor= %9.6f", i, Pars.detpolang[i], Pars.detdopfac[i]);
      printf ("1/%9.6f\n", 1.0 / Pars.detdopfac[i]);
    };

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

  GTID idbuf;
  sprintf (idbuf.IDSTR1, "GRETINA");
  sprintf (idbuf.IDSTR2, "GRETINA");
  idbuf.len = IDLEN;

  ctkROOTSetup ();

  /* allocate space for the crystal hits structure */

  track.payload = (PAYLOAD *) calloc (MAXTRACK, sizeof (PAYLOAD));
  track.gd = (GEBDATA *) calloc (MAXTRACK, sizeof (GEBDATA));

  /*--------------------------------*/
  /* loop and read the tracked data */
  /*--------------------------------*/

  printf ("\n");
  printf ("start reading data...\n");
  printf ("\n");

  st = getTrackedEvent (&track);

  while (st == 0 && Pars.evno < 2000000000)
    {

      if (Pars.evno % 10000 == 0)
        {
          fprintf (stdout, "track # %10i\n", Pars.evno);
          fflush (stdout);
        };

      if (Pars.nprint > 0)
        {

          sprintf (str, "trackedEvent_%3.3i.txt", Pars.evno);
          fp = fopen (str, "w");
          printEvent (fp, Pars.evno, &track);
          printEvent (stdout, Pars.evno, &track);
          fclose (fp);
        };
      Pars.nprint--;

      /* find the detector numbers as we need */
      /* them for various reasons below */

      findDetectorNumbers (&track);

      /*------------------------------*/
      /* Dopper correct the energies  */
      /*------------------------------*/


      /* find the GT mean time based on first hit */


      nmeanBGS_time = 0;
      nmeanS800_time = 0;
      nmeanGT_time = 0;
      meanBGS_time = 0;
      meanS800_time = 0;
      meanGT_time = 0;
      ptinp = track.payload;
      ptgd = track.gd;
      for (k = 0; k < track.n; k++)
        {

          if (ptgd->type == GEB_TYPE_BGS)
            {
              meanBGS_time += ptgd->timestamp;
              nmeanBGS_time++;
            };

          if (ptgd->type == GEB_TYPE_GT_MOD29 || ptgd->type == GEB_TYPE_S800PHYSDATA)
            {
              meanS800_time += ptgd->timestamp;
              nmeanS800_time++;
            };

          if (ptgd->type == GEB_TYPE_TRACK)
            {
              meanGT_time += ptgd->timestamp;
              nmeanGT_time++;
            };


          /* next */

          ptinp++;
          ptgd++;
        };

      /* find external TAC value */

      if (meanGT_time > 0)
        {
          meanGT_time /= nmeanGT_time;
        };

      if (nmeanS800_time > 0)
        {
          meanS800_time /= nmeanS800_time;
          EXT_tac = (int) (meanS800_time - meanGT_time);
//            printf("EXT_tac=%i; ", EXT_tac);
//            printf("meanS800_time=%i; ",meanS800_time);
//            printf("meanGT_time=%i;\n",meanGT_time);
        };
      if (nmeanBGS_time > 0)
        {
          meanS800_time /= nmeanS800_time;
          EXT_tac = (int) (meanBGS_time - meanGT_time);
        };

//          if (EXT_tac < 0 || EXT_tac > 2047)
//            {
//          printf("ooops: ");
//          printf("meanGT_time = %lli; meanBGS_time = %lli, diff = %i\n", meanGT_time,meanBGS_time,EXT_tac);
//              EXT_tac = 0;
//              fflush (stdout);
//            };


      /* do the binning */

//          bgs.nEsi = 0;
//          if (ProcessExternal (&track, EXT_tac, &bgs, &ntac, tac))

//            printf("ctkana passes (track->n=%i):\n",track.n);
      ctkROOTBin (&track, polang, doppler_factor, &bgs, EXT_tac);


      /* get next events tracked data points */

      st = getTrackedEvent (&track);
      if (st != 0)
        printf ("getTrackedEvent returned %i at event no %i\n", st, Pars.evno);

    };



  printf ("sorted %i tracked events\n", Pars.evno);

  /* write spectra into root file */

  ctkROOTexit ();

  /* statistics */

  printf ("sizeof(PAYLOAD)=%i\n", sizeof (PAYLOAD));
  printf ("sizeof (CLUSTER_INTPTS)=%i\n", sizeof (CLUSTER_INTPTS));

  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];
  printf ("total= %i\n", i1);

  printf ("\n");

  /* detector hit statistics */

  printf ("\n");
  i1 = 0;
  i2 = 0;
  for (i = 0; i < MAXDETNO; i++)
    if (dethit[i] > 10)
      {
        printf ("det %3i hits: %10i\n", i, dethit[i]);
        i1++;
        i2 += dethit[i];
      };
  printf ("we saw %i crystals with hits\n", i1);
  printf ("we saw %i hits\n", i2);
  printf ("\n");

  mean_pol /= mean_n;
  mean_azi /= mean_n;
  printf ("mean polar/azi angles: %6.2f %6.2f\n", (float) mean_pol, (float) mean_azi);

  /* done */

  exit (0);
}