/* dmp_msu_bgo.c */

/* 95/07/01 tl; first try */
/* 95/07/02 tl; added most of the bgo stuff */
/* 95/07/03 tl; one more time at it, most 1D spectra there */
/* 95/07/05 tl; added 2D */

#include <stdio.h>
#include <string.h>
#include "bgo.h"
#include "msu.h"

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

main(argc, argv)
  int             argc;
  char           *argv[];
{

  /* declarations */

  static int      st = 0, nn = 0, i, np = 0, j;
  int             max_blocks, max_prints;
  unsigned short int ev[MAXEV], evlen;
  int             err[20];

  unsigned short int mdbgo;	/* del BGO multiplicity */
  unsigned short int mibgo;	/* iso BGO multiplicity */
  unsigned short int mpbgo;	/* prompt BGO multiplicity */
  unsigned short int ph;	/* prompt BGO sum energy */
  unsigned short int lge, ege[NGE], geid[NGE];	/* ge energy array */
  unsigned short int lbgo, ebgo[NBGO], bgoid[NBGO];	/* BGO energy */
  short int       tge[NGE];	/* ge time array */
  unsigned short int tbgo[NBGO];/* BGO time array */

  /* spectra related to BGO setup */

  char            str[100];
  float           rsp[M13];

  int             sege[NGE][M13];
  int             sum[M13];
  int             stge[NGE][M13];
  int             patge[NGE];

  int             sebgo[NBGO][M11];
  int             stbgo[NBGO][M11];
  int             patbgo[NBGO];

  int             m_del[NBGO], m_iso[NBGO], m_pro[NBGO];
  int             sumh[M11];

  int            *mat;
  int            *malloc();
  char           *p, *strstr();
  int             s2d, no_files, files_read = 0;
  FILE           *fp, *fopen();
  short int       e1, e2;

  /* help */

  if (argc != 5)
  {
    printf("use: dmp_msu <#blocks> <#events to print> <2D|no2D> <#files>\n");
    return (0);
  };

  /* hello */

  printf("dmp_msu.c running on ");
  time_stamp();

  /* get parameters */

  max_blocks = atoi(argv[1]);
  max_prints = atoi(argv[2]);

  if ((p = strstr(argv[3], "2D")) != NULL)
    s2d = 1;
  else if ((p = strstr(argv[3], "2d")) != NULL)
    s2d = 1;
  else
    s2d = 0;
  if ((p = strstr(argv[3], "no")) != NULL)
    s2d = 0;
  if ((p = strstr(argv[3], "NO")) != NULL)
    s2d = 0;


  no_files = atoi(argv[4]);


  /* zap spectra */

  for (i = 0; i < 20; i++)
    err[i] = 0;
  for (i = 0; i < M13; i++)
    for (j = 0; j < NGE; j++)
    {
      sege[j][i] = 0;
      stge[j][i] = 0;
    }

  /* allocate the 2D matrix if user asked for it */

  if (s2d)
  {
    mat = malloc(MATLEN * MATLEN * sizeof(int));
    if (mat == NULL)
    {
      printf("not enough memory for matrix\n");
      exit(0);
    } else
    {
      printf("allocated %i bytes for matrix\n", MATLEN * MATLEN * sizeof(int));

      /* zap it just in case. should be zero but... */

      for (i = 0; i < MATLEN * MATLEN; i++)
	*(mat + i) = 0;

    }
  };

  /* sort */

  printf("sorting started @ ");
  time_stamp();

  while (st >= 0 && err[1] <= max_blocks && files_read < no_files)
  {

    /* get an event and count */

    st=bgo_ev(&evlen, ev,
	     &mdbgo, &mibgo, &mpbgo,
	     &ph,
	     &lge, ege, geid, tge,
	     &lbgo, ebgo, bgoid, tbgo);

    nn++;

    /* count files */

    if (st == FILEMARK)
      files_read++;


    /* keep log of the errors/counts */

    if ((st > -1) && (st < 20))
      err[st]++;
 
    /* process good event */

    if (st == 0)
    {

      /* print some events */

      if (np < max_prints)
	{
	  np++;
	  printf("event print @ call # %i:\n", nn);

	  /* raw format via routine in msu.c */

	  print_event(&evlen, ev);
          fflush(stdout);

	  /* interpreted BGO data */

	  print_bgo(&evlen, ev,
		    &mdbgo, &mibgo, &mpbgo,
		    &ph,
		    &lge, ege, geid, tge,
		    &lbgo, ebgo, bgoid, tbgo);
          fflush(stdout);

	};

      /* bin germanium energies and times */

      for (i = 0; i < lge; i++)
      {
	sege[geid[i]][ege[i]]++;
	stge[geid[i]][tge[i]]++;
      }

      /* bin BGO energies and times */

      for (i = 0; i < lbgo; i++)
      {
	sebgo[bgoid[i]][ebgo[i]]++;
	stbgo[bgoid[i]][tbgo[i]]++;
      }

      /* bin BGO multiplicities and sum energy */

      m_del[mdbgo]++;		/* delayed */
      m_iso[mibgo]++;		/* isomeric */
      m_pro[mpbgo]++;		/* prompt */
      sumh[ph]++;		/* prompt summed energy in BGO ball */

      /* 2D */

      if (s2d)
	if (lge >= 2)
	{
	  for (i = 0; i < lge; i++)
	    for (j = i + 1; j < lge; j++)
	    {
	      e1 = ege[i] / 2;
	      if (e1 < MATLEN)
	      {
		e2 = ege[j] / 2;
		if (e2 < MATLEN)
		{

		  MAT_INCREMENT(e1, e2, 1);
		  MAT_INCREMENT(e2, e1, 1);
		}
	      }
	    }


	}
    }
  }

  printf("\ndone sorting @ ", st);
  time_stamp();

  err_print(st);
  printf("vas the last error code! (%i)\n\n", st);

  /* statistics */

  printf("counts and errors encountered:\n");
  for (i = 0; i < 20; i++)
    if (err[i] > 0)
    {
      err_print(i);
      printf(" %10i\n", err[i]);
    };

  if (nn>10)
  {

  /* write ge related spectra out */

  for (i = 0; i < M13; i++)
    sum[i] = 0;
  for (i = 0; i < NGE; i++)
    patge[i] = 0;
  for (j = 1; j < NGE; j++)
  {
    for (i = 0; i < M13; i++)
    {
      rsp[i] = (float) sege[j][i];
      patge[j] += sege[j][i];
      sum[i] += sege[j][i];
    }
    i = M13;
    sprintf(str, "ege_%2.2i.spe", j);
    wr_spe(str, &i, rsp);
    printf("[%s]\n", str);

  }
  for (i = 0; i < M13; i++)
    rsp[i] = (float) sum[i];
  i = M13;
  sprintf(str, "ege_sum.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);
  rsp[0] = 0;
  for (i = 1; i < NGE; i++)
    rsp[i] = (float) patge[i];
  i = NGE;
  sprintf(str, "hitge.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);
  for (j = 1; j < NGE; j++)
  {
    for (i = 0; i < M13; i++)
      rsp[i] = (float) stge[j][i];
    i = M13;
    sprintf(str, "tge_%2.2i.spe", j);
    wr_spe(str, &i, rsp);
    printf("[%s]\n", str);

  }

  /* write bgo related spectra out */

  for (i = 0; i < M11; i++)
    sum[i] = 0;
  for (i = 0; i < NBGO; i++)
    patbgo[i] = 0;
  for (j = 1; j < NBGO; j++)
  {
    for (i = 0; i < M11; i++)
    {
      rsp[i] = (float) sebgo[j][i];
      patbgo[j] += sebgo[j][i];
      sum[i] += sebgo[j][i];
    }
    i = M11;
    sprintf(str, "ebgo_%2.2i.spe", j);
    wr_spe(str, &i, rsp);
    printf("[%s]\n", str);

  }
  for (i = 0; i < M11; i++)
    rsp[i] = (float) sum[i];
  i = M11;
  sprintf(str, "ebgo_sum.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);
  rsp[0] = 0;
  for (i = 1; i < NBGO; i++)
    rsp[i] = (float) patbgo[i];
  i = NBGO;
  sprintf(str, "hitbgo.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);
  for (j = 1; j < NBGO; j++)
  {
    for (i = 0; i < M11; i++)
      rsp[i] = (float) stbgo[j][i];
    i = M11;
    sprintf(str, "tbgo_%2.2i.spe", j);
    wr_spe(str, &i, rsp);
    printf("[%s]\n", str);

  }

  /* ball multiplicity and summed prompt energy */

  rsp[0] = 0;

  for (j = 1; j < NBGO; j++)
    rsp[j] = (float) m_del[j];
  i = M11;
  sprintf(str, "m_del.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);

  for (j = 1; j < NBGO; j++)
    rsp[j] = (float) m_iso[j];
  i = M11;
  sprintf(str, "m_iso.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);

  for (j = 1; j < NBGO; j++)
    rsp[j] = (float) m_pro[j];
  i = M11;
  sprintf(str, "m_pro.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);

  for (j = 1; j < M11; j++)
    rsp[j] = (float) sumh[j];
  i = M11;
  sprintf(str, "sumh.spe");
  wr_spe(str, &i, rsp);
  printf("[%s]\n", str);

  /* write 2D out */

  if (s2d)
  {
    fp = fopen("bgo.mpg", "w");
    st = fwrite(mat, MATLEN * MATLEN * sizeof(int), 1, fp);
    fclose(fp);
    if (st == 1)
      printf("2D: %i bytes --> 'bgo.mpg'\n", MATLEN * MATLEN * sizeof(int));
    else
      printf("failed to write matrix!\n");
  };

  /* done */

  };

  printf("** done **\n");
  return (0);
}

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