/* $Id: GSSort.cxx,v 1.330 2011/03/07 20:09:10 tl Exp $ */

/* GS data sorting program for ROOT */
/* the meanest data sifter on the block! */
/* by Torben Lauritsen, ANL */

/* set SIXTYFOUR to 1 if you are on a 64 bit */
/* machine and 0 if on a 32 bit machine */

#define SIXTYFOUR 0

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#ifdef SOLARIS
#include <strings.h>
#endif
#ifdef LINUX
#include <string.h>
#endif
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <signal.h>
#include <time.h>
#include <stddef.h>

#include "Rtypes.h"
#include "TROOT.h"
#include "TFile.h"
#include "TRandom.h"
#include "TH1.h"
#include "TH2.h"
#include "TObjArray.h"
#include "TObject.h"
#include "TKey.h"
#include "TSystem.h"
#include "TCutG.h"
#include "TTree.h"

#include "TMapFile.h"

#include "gsII.h"
#include "GSSort.h"
#include "GSudpReceiver.h"

/* Don's tree stuff */

#include "UserInclude.h"

/* global sort parameters with defaults */


#include "GSSort_global_declares.h"

unsigned int    bgoseghit[NGE][7];
unsigned int    gehits[NGE];
int             GSudpPort = 1101;
unsigned int    nEvents = 1000000000;
int             firstEvent = 0;
char            InputSource[STRLEN] = "/dev/rmt/1mbn";
char            ROOTFile[STRLEN] = "GSSort.root";
char            treeFile[STRLEN] = "GSFMATree.root";
int             UseRootFile = FALSE, SizeShareMemFile = FALSE;
char            ShareMemFile[STRLEN] = "GSSort.map";
int             UseShareMemFile = FALSE;
char            ROOTFileOption[STRLEN] = "UNDEFINED";
int             NumToPrint = 25;
float           Beta = 0;
int             InputSrc = NOTDEF;
int             FirstFile = NOTDEF;
int             Nfiles = 100;
int             WriteRawFera = FALSE;
int             NFERAErrToPrint = 10;
int             mpevOverride = FALSE;
int             DoNotInterpExt = FALSE;

int             PrintRawPEv = FALSE;

int             DumpEvery = 10;
Int_t           ComPressLevel = NOTDEF;
#ifdef SIXTYFOUR
 unsigned long long int    StartMapAddress = 0;
#else
 unsigned int    StartMapAddress = 0;
#endif
int             CountPevLevCheck = 0;
int             UpdateRootFile = 0;
int             sub = FALSE;

/* other globals */

int             WeWereSignalled = FALSE;	/* signal  */
int             tStart;
char            PevsInitialized = FALSE;
int             NprintEvNo = 0;
int             InputSourceOnCommandLine = 0;
int             HaveRootFileName = 0;

double          angle[NGE], DopCorFacHires[NGE], DopCorFac[NGE], ACFac[NGE];

int             set_gsII_indrive(char *);
int             SetupUdpReceiver(int);
void            CheckNoArgs(int, int, char *);
char            CommandFileName[STRLEN] = "GSSort.command";

/* macros */

#include "GSGTMacros.h"

/* allow for user functions here */

#include "UserFunctions.h"

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

void
printBGOSegmentHits(void)
{

  /* declarations */

  int             i, j;
  double          sum = 0, val = 0;

  /* headline */

  printf("\n");
  printf("BGO segment hit statistics:\n");
  printf("\n");

  /* go list */

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

    /* first sum up all 7 channels */

    sum = 0;
    for (j = 0; j < 7; j++)
      sum += (double) bgoseghit[i][j];

    printf("BGO%3.3i: %9.0f -> ", i, (float) sum);

    /* then the segments in percent */

    if (sum > 1.0)
    {

      for (j = 0; j < 7; j++)
      {
	val = 100.0 * (double) bgoseghit[i][j] / sum;
	printf("%3i%% ", (int) (val + 0.5));
      };

      /* bgo to germanium hits */

      if (gehits[i] > 1)
      {
	val = 100.0 * (double) bgoseghit[i][j] / (double) gehits[i];
	printf("; Ge/BGO=%3i%% ", (int) (val + 0.5));
      } else
	printf("; no ge counts!?");

    } else
    {
      printf("-- BGO counts=%i, Ge counts=%i", (int)sum, (int)gehits[i]);
    };  


    /* done with this detector */

    printf("\n");

  };

  printf("\n");
  printf("__ even counts would yield ~14%% in each segment\n");
  printf("\n");

};


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

void
signal_catcher(int sigval)
{
  int             time_stamp();
  printf("\nGSSort/GSacq received signal <%i> on ", sigval);
  time_stamp();
  WeWereSignalled = TRUE;
  fflush(stdout);
}

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

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

  /*--------------*/
  /* declarations */
  /*--------------*/

  int             j, i, HaveChatFile = 0;
  char           *p;
  char            ChatFileName[STRLEN];
  int             GSacq(char *);
  int             time_stamp();
  char            str2[STRLEN], str3[STRLEN], str4[STRLEN];

  /* inform of version of root used */

  printf("from %s on %s\n", __FILE__, __DATE__);
  printf("RCS version $Id: GSSort.cxx,v 1.330 2011/03/07 20:09:10 tl Exp $\n");

  /* delay a few sec (wait for logging window to come up) */

  printf("VVV sleep 5 sec (wait for logger window)...\n");
  fflush(stdout);
  sleep(5);

  /*------*/
  /* help */
  /*------*/

  if (argc < 2)
  {
    printf("\n");
    printf("use: %s -chat file [-help] [-version] [-input disk|net|tape file|drive|port UPDATE|RECREATE  -rootfile name.root]\n", argv[0]);
    printf("\n");
    return (0);
  };

  /*--------------------*/
  /* parse command line */
  /* and call GSacq     */
  /*--------------------*/

  j = 1;			/* current command line arg position */

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

      if ((p = strstr(argv[j], "-version")) != NULL)
      {
	printf("$Id: GSSort.cxx,v 1.330 2011/03/07 20:09:10 tl Exp $\n");
	printf("This version has Float 2D matrices\n");
	exit(0);
      } else if ((p = strstr(argv[j], "-help")) != NULL)
      {
	printf("\n");
	printf("GSSort is documented on the WWW, URL:\n");
	printf("\n");
	printf("  http://www.phy.anl.gov/gs/doc/GSSort\n");
	printf("\n");
	exit(0);

      } else if ((p = strstr(argv[j], "-input")) != NULL)
      {

	/* check that user specified enough arguments */

	j++;
	if ((argc - j) < 3)
	{
	  printf("you must specify -input disk|net|tape file|drive|port UPDATE|RECREATE\n");
	  exit(0);
	}
	/* input source will be specified on command line. */
	/* Set flag to prevent processing of */
	/* input source in chat script  */

	InputSourceOnCommandLine = 1;

	/* determine input source */

	strcpy(str2, argv[j++]);
	strcpy(str3, argv[j++]);
	strcpy(ROOTFileOption, argv[j++]);

	printf("\n");
	printf("*** GSSort will take \"input source\" from command line\n");

	if (strcmp("tape", str2) == 0)
	{
	  printf("will take input from tape\n");
	  set_gsII_indrive(str3);
	  InputSrc = TAPE;
	  fflush(stdout);
	} else if (strcmp("disk", str2) == 0)
	{
	  printf("will take input from disk\n");
	  set_gsII_indrive(str3);
	  InputSrc = DISK;
	  fflush(stdout);
	} else if (strcmp("net", str2) == 0)
	{
	  printf("will take input from UDP sender\n");
	  GSudpPort = atoi(str3);
	  SetupUdpReceiver(GSudpPort);
	  printf("__using port %i\n", GSudpPort);
	  InputSrc = NET;
	  fflush(stdout);
	} else
	{
	  printf("unknown input option: %s\n", str2);
	  printf("aborting\n");
	  fflush(stdout);
	  exit(0);
	};

	if ((p = strstr(ROOTFileOption, "UPDATE")) != NULL)
	{
	  UpdateRootFile = TRUE;
	  printf("will update root file\n");
	} else if ((p = strstr(ROOTFileOption, "RECREATE")) != NULL)
	{
	  UpdateRootFile = FALSE;
	  printf("will recreate root file\n");
	} else
	{
	  printf("option <%s> not valid!\n", str4);
	  printf("valid options are: RECREATE or UPDATE\n\n");
	  exit(0);
	};

	printf("\n");

      } else if ((p = strstr(argv[j], "-chat")) != NULL)
      {
	j++;
	strcpy(ChatFileName, argv[j++]);
	printf("will read sorting instructions from chatfile: %s\n", ChatFileName);
	system("ulimit -a");
	HaveChatFile = 1;
      } else if ((p = strstr(argv[j], "-rootfile")) != NULL)
      {
	j++;
	strcpy(ROOTFile, argv[j++]);
	printf("rootfile name specified on command line\n");
	printf("__will store spectra in rootfile: %s\n", ROOTFile);
	HaveRootFileName = 1;
	UseRootFile = 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);
	}
	printf("]\non ");
	time_stamp();
	exit(0);

      }
    };

  /* now start the sorter */

  if (HaveChatFile)
    GSacq(ChatFileName);
  else
  {
    printf("you must specify a chat script\n");
    exit(0);
  }



};


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

int
GSacq(char *ChatFileName)
{

  /*--------------*/
  /* declarations */
  /*--------------*/

  /* root spectra etc */

#include "GSSort_root_spectra.h"


  TH1D           *PEv1D[NFERA];
  TH2F           *PEv2D[NFERA];
  TH2F           *GamPEv2D[NFERA];
  TH2F           *sDevMatrix[MAXDEVMATRIX];
  TH2F           *sHit2DMat[MAXDEVMATRIX];
  TH2F           *GeTEmat[NGE];

  TH1D           *evrate, *ehirate;
  TH1D           *rate_sp[NGE];
  TH1D           *datarate;

  TH2F           *TMul_ehi;
  TH2F           *ehi_ehi[MAX2DS];

  TH1            *htemp;

  TList          *wlist;
  TList          *zlist;
  TIterator      *hiterator;

  TMapFile       *mfile, *m;
  TFile          *f1;

  TH1D           *tmpTH1D = NULL;
  TH2F           *tmpTH2F = NULL;


  /* misc */

  static int      bincode[16] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768};
  char            str1[STRLEN], str2[STRLEN];
  int             i, j, k, l, st = 0, nn, nret;
  int             ok;
  int             i0, i1, i2, i3;
  float           r1;
  FILE           *fp;
  char           *pc;
  unsigned int    seed = 0;
  int             exa;
  unsigned int    CurEvNo = 0;
  int             CurFile = 0;
  int             eov, oldst;
  int             lenclean, lendirty, lentotal;
  int             MultOK, nMultOK = 0;
  time_t          t1, t2, tstart, tnow, tdmplast, tdmp = 0;
  double          EvTime = 0;
  char            str[STRLEN], spname[STRLEN];
  char           *p;
  int             IsResidue = 0, nres = 0, nppac = 0, cpev;
  unsigned int    rnseed;
  int             FERAPEvlookup[MAXVSN][MAXFERACH];
  double          d1, d2, d3;
  double          emax, rn;
  int             okge, okbgo, ok1, ok2;
  unsigned short int tmp_lc, tmp_ld, tmp_lb;

  int             nxx = 0;
  float           r2;
  unsigned short int start_ttL , start_ttM, start_ttH;
  int             firstevent=1;
  double          runTime;

  /*---------------------*/
  /* function prototypes */
  /*---------------------*/

  int             time_stamp();
  int             get_gsII_ev(struct EVHDR *, struct EVENT *, unsigned int *);
  int             get_disk_ev(struct EVHDR *, struct EVENT *, unsigned int *);
  int             get_net_ev(struct EVHDR *, struct EVENT *, unsigned int *);
  int             GSSort_read_chat(char *);
  int             print_event(struct EVHDR, struct EVENT, int, int, int, float, float, int);
  int             decode_fera(unsigned int *, unsigned int *,
			                      unsigned int *, unsigned int *, unsigned int *, unsigned int *,
			                      double *, double *);
  int             decode_isomer(unsigned int *, unsigned int *,
				                unsigned int *, unsigned int *, unsigned int *, unsigned int *,
				                unsigned short int);
  int             fera_print(int, unsigned int *, unsigned int *,
			                     unsigned int *, unsigned int *);
  int             get_a_seed(unsigned int *);
  int             status_tape(int);
  int             rew_tape(int);
  int             fskp_tape(int, int);
  int             CloseUdpReceiver();
  void            WriteStatistics(unsigned int, int, int *, double);
  TCutG          *rd2dwin(Char_t *);
  void            ApplyAllConditions(void);
  void            FindConditions(void);
  void            FindParentPEv(int);
  void            SetFERAErrorPrint(int);
  void            PrintFERATypes();
  void            CountConditions();
  void            mpevOP_reg();
  void            mpevOP_special();
  void            MarkCountingPevs();
  TH1D           *mkTH1D(char *, char *, int, double, double);
  TH2F           *mkTH2F(char *, char *, int, double, double, int, double, double);
  void            SetBeta();
  void            ReportBufPosErrors(void);
  void            PrintCondToApply();
  void            wrIsomerStatistics(int);
  int             create_isomerdata(struct EVHDR *, struct EVENT *, unsigned int *);
  void            DoPrimaryGating(void);
  void            ResetGSII_bufpCounters();
  int             nxx1 = 200;
  void            ZapFERAErrors(void);

  /* allow user to declare variables here */

  printf("\n");
  printf("executing UserDeclare.h code\n");
  printf("\n");
#include "UserDeclare.h"

  /*-------*/
  /* Hello */
  /*-------*/

  printf("\n");
  printf("GSacq running on: ");
  time_stamp();
  printf("$Id: GSSort.cxx,v 1.330 2011/03/07 20:09:10 tl Exp $\n");
  printf("\n");

  /*------------------------------------------*/
  /* name the Gammasphere entries of the pev  */
  /* so they can be recognized in chat script */
  /*------------------------------------------*/

  printf("\n");
  printf("The following are standard GS pevs\n");
  printf("__which will be available in this sort\n");
  printf("\n");
  fflush(stdout);

  sprintf(PEv[GAMMASPHERE].name, "gammasphere");
  printf("%s\n", PEv[GAMMASPHERE].name);
  sprintf(PEv[LEN_CLEAN].name, "mclean");
  printf("%s\n", PEv[LEN_CLEAN].name);
  sprintf(PEv[LEN_DIRTY].name, "mdirty");
  printf("%s\n", PEv[LEN_DIRTY].name);
  sprintf(PEv[LEN_BGO].name, "mbgo");
  printf("%s\n", PEv[LEN_BGO].name);
  sprintf(PEv[LENCLEAN].name, "lenclean");
  printf("%s\n", PEv[LENCLEAN].name);
  sprintf(PEv[LENDIRTY].name, "lendirty");
  printf("%s\n", PEv[LENDIRTY].name);
  sprintf(PEv[LENTOTAL].name, "lentotal");
  printf("%s\n", PEv[LENTOTAL].name);
  sprintf(PEv[SUMGE].name, "sumge");
  printf("%s\n", PEv[SUMGE].name);
  sprintf(PEv[SUMBGO].name, "sumbgo");
  printf("%s\n", PEv[SUMBGO].name);
  sprintf(PEv[CLOCK1].name, "clock1");
  printf("%s\n", PEv[CLOCK1].name);
  sprintf(PEv[CLOCK2].name, "clock2");
  printf("%s\n", PEv[CLOCK2].name);
  sprintf(PEv[TAC2].name, "tac2");
  printf("%s\n", PEv[TAC2].name);
  sprintf(PEv[LENEXT].name, "lenext");
  printf("%s\n", PEv[LENEXT].name);
  fflush(stdout);

  for (i = 0; i < NGSPEV; i++)
  {
    sprintf(PEv[i + GSPEVL1].name, "egehi%2.2i", i);
    sprintf(PEv[i + GSPEVL2].name, "tgerf%2.2i", i);
  };
  printf("\n");
  fflush(stdout);

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

  printf("\n");
  printf("initializing\n");
  printf("\n");
  fflush(stdout);

  /* time offsets */

  for (i = 0; i < NGE; i++)
  {
    tGEoffset[i] = 0;
    tBGOoffset[i] = 0;
  };

  printf("done initalizing\n");
  fflush(stdout);

  /*------------------*/
  /* read chat script */
  /*------------------*/

  GSSort_read_chat(ChatFileName);
  PrintFERATypes();

  /* make sure we mark all counting pevs */

  MarkCountingPevs();

  /*---------------*/
  /* sanety checks */
  /*---------------*/

  if (InputSrc == NOTDEF)
  {
    printf("you must specify an input source!\n");
    printf("quitting...\n");
    exit(1);
  };

  if (UseShareMemFile && UseRootFile)
  {
    printf("you cannot use shared memory and a root disk\n");
    printf("at the same time!\n");
    exit(1);
  };

  /* force user to declare intension with root file */
  /* so I can't be blamed for any overwrites!!      */

  if (!UseShareMemFile)
    if ((p = strstr(ROOTFileOption, "UNDEFINED")) != NULL)
    {
      printf("for root files you must specify either:\n");
      printf("\n");
      printf("    rootfileoption RECREATE\n");
      printf("or\n");
      printf("    rootfileoption UPDATE\n");
      printf("\n");
      printf("please modify your chat file and try again\n");
      exit(-1);
    };

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

  hdr.len = 0;			/* force record read first time */

  get_a_seed(&seed);		/* initialize random number generator */
  srand(seed);

  for (i = 0; i < NERR; i++)
    GSSortError[i] = 0;

  CurEvNo = 0;
  NprintEvNo = 0;

  if (!UseShareMemFile)
  {
    DumpEvery = 1000000000;
    printf("\n");
    printf("_since rootfile: setting `DumpEvery` to infinity..!\n");
    printf("\n");
  };

  /* delete any command file */

  sprintf(str, "\\rm -f %s", CommandFileName);
  system(str);
  printf("%s\n", str);

  /*--------------------------------------*/
  /* read old rootfile ? - continue sort? */
  /*--------------------------------------*/

  if (UpdateRootFile)
  {

    /* check here whether the old root file exists */

    fp = fopen(ROOTFile, "r");
    if (fp == NULL)
    {
      printf("could not open old rootfile: %s\n", ROOTFile);
      printf("the old rootfile must exist if you  \n");
      printf("want to use the UPDATE option\n");
      printf("aborting...\n");
      exit(0);
    };
    fclose(fp);

    /* read in old root file */

    f1 = NULL;
    f1 = new TFile(ROOTFile, "UPDATE");
    printf("read old root file <%s>\n", ROOTFile);
    if (!f1->IsOpen())
    {
      printf("could not open file....\n\n");
      exit(-1);
    };
    printf("base=<%s>\n", f1->GetPath());
    f1->Print();

  };

  /*---------------------------------------*/
  /* list pseudo event vector assignments */
  /*---------------------------------------*/

  printf("\n");
  printf("GAMMASPHERE pseudoevent vector assignments:\n");

  for (i = 0; i < GSPEVL1; i++)
    printf("pev %3.3i is <%s>\n",
	   i, PEv[i].name);
  printf("\n");

  for (k = 0; k < MAXVSN; k++)
    for (l = 0; l < MAXFERACH; l++)
      FERAPEvlookup[k][l] = -1;

  printf("FERA pseudoevent vector assignments:\n");
  for (i = NPEVMIN; i < LenPEv; i++)
    if (PEv[i].vsn != 0)
    {
      printf("pev  name: <%s>[#%3.3i] is FERA VSN: %3.3i/CH:%2.3i; ",
	     PEv[i].name, i, PEv[i].vsn, PEv[i].ch);
      printf("tresh:%4i ", PEv[i].thresh);
      if (PEv[i].calib)
	printf("calib: %f +%f*ch\n", (float) PEv[i].offset, (float) PEv[i].gain);
      else
	printf("no calibration\n");

      /* check */

      if (FERAPEvlookup[PEv[i].vsn][PEv[i].ch] == -1)
	FERAPEvlookup[PEv[i].vsn][PEv[i].ch] = i;
      else
      {
	printf("ooops\n");
	exit(-1);
      };

    };
  printf("\n");

  /* check lookup table */

  for (i = 0; i < MAXVSN; i++)
    for (j = 0; j < MAXFERACH; j++)
      if (FERAPEvlookup[i][j] > 0)
	if (PEv[FERAPEvlookup[i][j]].vsn != i || PEv[FERAPEvlookup[i][j]].ch != j)
	{
	  printf("inconsistent: i=%i,j=%i\n", i, j);
	  printf("FERAPEvlookup[i][j]=%i\n", FERAPEvlookup[i][j]);
	  printf("PEv[FERAPEvlookup[i][j]].vsn=%i ", PEv[FERAPEvlookup[i][j]].vsn);
	  printf("PEv[FERAPEvlookup[i][j]].ch=%i\n", PEv[FERAPEvlookup[i][j]].ch);
	  exit(-1);
	};

  printf("other pseudoevent vector assignments:\n");
  for (i = NPEVMIN; i < LenPEv; i++)
    if (PEv[i].vsn == 0)
    {
      printf("pev %3.3i, name: <%s>", i, PEv[i].name);
      if (PEv[i].counting)
	printf(" [COUNTER PEV]\n");
      else
	printf("\n");
    };
  printf("\n");
  printf("\n");

  if (NPEvCond > 0)
  {

    printf("\n");
    printf("Conditions->\n");
    printf("\n");
    printf("the length of the condition vector is %i\n", NPEvCond);

    /*--------------------*/
    /* list 1D conditions */
    /*--------------------*/

    printf("\n");
    printf("1D conditions\n");
    printf("\n");
    i1 = 0;
    for (i = 0; i < NPEvCond; i++)
      if (PEvCond[i].filled)
	if (PEvCond[i].type == 1)
	{
	  printf("1D cond: %3i: on pev %3i, <%s>, ", i, PEvCond[i].pev1, PEvCond[i].name);
	  printf("[%s]\n", PEvCond[i].strrange);
	  i1++;
	};
    printf("\n");
    printf("__there are %i 1D conditions\n", i1);

    /*--------------------*/
    /* list 2D conditions */
    /*--------------------*/

    printf("\n");
    printf("2D conditions\n");
    printf("\n");
    i1 = 0;
    for (i = 0; i < NPEvCond; i++)
      if (PEvCond[i].filled)
	if (PEvCond[i].type == 2)
	{
	  printf("2D cond: %3i: on pev %3i and %3i, ", i, PEvCond[i].pev1, PEvCond[i].pev2);
	  printf("<%s>,  with 2Dwinname: %s\n", PEvCond[i].name, PEvCond[i].fn);
	  /* PEvCond[i].win->Print(); */
	  i1++;
	};
    printf("\n");
    printf("__there are %i 2D conditions\n", i1);
    printf("\n");

    /*-------------------------*/
    /* find all derived pevs   */
    /* and apply the necessary */
    /* inherited conditions    */
    /*-------------------------*/

    printf("\n");
    for (i = NPEVMIN; i < LenPEv; i++)
      if (PEv[i].vsn == 0)
      {
	printf("derived pev: [#%i,name:%s]\n", i, PEv[i].name);

	/* find parent pev */

	FindParentPEv(i);

      };
    printf("\n");

    /*-----------------------------------*/
    /* list ALL conditions to be applied */
    /* and mark inherited conditions     */
    /*-----------------------------------*/

    PrintCondToApply();

  };				/* if (NPEvCond > 0) */

  /*---------------------*/
  /* list all time masks */
  /*---------------------*/

  printf("\n");
  printf("time masks->\n");
  for (i = 0; i < NTimeMasks; i++)
    printf("%3i: time masks <%s> has mask {%s}\n", i, TMask[i].name, TMask[i].strrange);
  printf("\n");

  /*----------------------*/
  /* set up shared memory */
  /*----------------------*/

  if (UseShareMemFile)
  {
    printf("\n");

    if (StartMapAddress != 0)
    {
    
      printf("StartMapAddress will be set manually\n");
           
    } else
    {
    
      /* automatically find address (suggestion by Stefanos) */
      /* Note: this does not seem to work on SUNs... */
      
      printf("Auto find StartMapAddress...\n");

      /* make dummy map and find StartMapAddress */

      m = TMapFile::Create("dummy.map", "recreate", SizeShareMemFile);
#ifdef SIXTYFOUR
      StartMapAddress =  (unsigned long long int) m->GetMmallocDesc();  
#else
      StartMapAddress =  (unsigned int) m->GetMmallocDesc();
#endif
      printf("shared mem start address determined to be 0x%8.8x\n", StartMapAddress);
      m->Print();

      /* close and remove dummy map file */

      m->Close();
      gSystem->Exec("\\rm dummy.map");      
    
    };
        
    /* set the StartMapAddress address */
    
    TMapFile::SetMapAddress((Long_t) StartMapAddress);
    printf("shared mem start address set to 0x%8.8x\n", StartMapAddress);

    mfile = TMapFile::Create(ShareMemFile, "RECREATE", SizeShareMemFile, "GS shared mem");
    if (mfile == NULL)
    {
      printf("failed to create TMapFile\n");
      exit(-1);
    };
    
    printf("shared memory [%s] created, size: %i bytes\n", ShareMemFile, SizeShareMemFile);
    fflush(stdout);
    mfile->Print();
    printf("\n");

  };

  /*---------------*/
  /* position tape */
  /*---------------*/

  if (InputSrc == TAPE)
    if (FirstFile >= 0)
    {
      printf("\n");
      printf("positioning tape... current:\n");
      fflush(stdout);
      if ((exa = open(InputSource, O_RDONLY, 0)) == -1)
      {
	printf("could not opentape drive\n");
	exit(1);
      };
      status_tape(exa);
      fflush(stdout);
      rew_tape(exa);
      if (FirstFile > 0)
	fskp_tape(exa, FirstFile);
      printf("new position:\n");
      fflush(stdout);
      status_tape(exa);
      close(exa);
      fflush(stdout);
    };

  /*----------------------------------------*/
  /* get GS angles and find doppler factors */
  /*----------------------------------------*/

  /* get the default angles */

#include "GSII_angles.h"
  printf("default GS angles loaded from [GSII_angles.h]\n");

  /* optionally read in users angles */

  if (HaveAngleFileName)
  {

    /* attempt to open [optional] angle file */

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

    /* read content and set angles */

    pc = fgets(str1, STRLEN, fp);
    while (pc != NULL)
    {

      sscanf(str1, "%i %f", &i1, &r1);
      angle[i1] = (double) r1;
      printf("user specified: det %3.3i at angle %9.2f\n", i1, (float) angle[i1]);
      fflush(stdout);

      /* read next line */

      nn++;
      pc = fgets(str1, STRLEN, fp);

    };
    printf("read %i angles from %s\n", nn, AngleFileName);

    fclose(fp);

  }
  /* calculate beta=v/c related variables */

  SetBeta();

  /*---------------------------*/
  /* ring detectors angles etc */
  /*---------------------------*/

  ring_angle[0] = 0;
  ring_angle[1] = 17.27;
  ring_angle[2] = 31.72;
  ring_angle[3] = 37.38;
  ring_angle[4] = 50.07;
  ring_angle[5] = 58.28;
  ring_angle[6] = 69.82;
  ring_angle[7] = 79.19;
  ring_angle[8] = 80.71;
  ring_angle[9] = 90.0;
  ring_angle[10] = 99.29;
  ring_angle[11] = 100.81;
  ring_angle[12] = 110.18;
  ring_angle[13] = 121.78;
  ring_angle[14] = 129.93;
  ring_angle[15] = 142.62;
  ring_angle[16] = 148.28;
  ring_angle[17] = 162.73;

  for (i = 0; i < NGE; i++)
    angno[i] = 0;
  for (i = 1; i < NGE; i++)
    switch ((int) (angle[i] + 0.5))
    {
    case 17:
      angno[i] = 1;
      break;
    case 32:
      angno[i] = 2;
      break;
    case 37:
      angno[i] = 3;
      break;
    case 50:
      angno[i] = 4;
      break;
    case 58:
      angno[i] = 5;
      break;
    case 70:
      angno[i] = 6;
      break;
    case 79:
      angno[i] = 7;
      break;
    case 81:
      angno[i] = 8;
      break;
    case 90:
      angno[i] = 9;
      break;
    case 99:
      angno[i] = 10;
      break;
    case 101:
      angno[i] = 11;
      break;
    case 110:
      angno[i] = 12;
      break;
    case 122:
      angno[i] = 13;
      break;
    case 130:
      angno[i] = 14;
      break;
    case 143:
      angno[i] = 15;
      break;
    case 148:
      angno[i] = 16;
      break;
    case 163:
      angno[i] = 17;
      break;
    default:
      printf("angle %f not found\n", angle[i]);
      break;
    };
  printf("polar angles collected\n");
  fflush(stdout);

  /* print what detectors are in what angle */

  printf("\ndetector map:\n");
  for (i = 1; i <= 17; i++)
  {
    k = 0;
    printf("ang # %2i; ", i);
    fflush(stdout);
    for (j = 0; j < NGE; j++)
      if (angno[j] == i)
      {
	if (k == 0)
	{
	  k = 1;
	  printf("{%3i deg}> ", (int) (angle[j] + 0.5));
	}
	printf("[%3i]", j);
	fflush(stdout);
      }
    printf("\n");
  };
  printf("\n");

  /*-----------------------------------*/
  /* make root histograms if we didnt  */
  /* read in an old root file --------- */
  /*-----------------------------------*/

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

    /* hi res spectra */

    sprintf(str1, "ehi%3.3i", i);
    sprintf(str2, "det %3.3i hi res @ angle %6.2f", i, (float) angle[i]);
    ehi[i] = mkTH1D(str1, str2, L14BITS, 1, M14BITS);

    /* eside spectra */

    sprintf(str1, "eside%3.3i", i);
    sprintf(str2, "det %3.3i side channel @ angle %6.2f", i, (float) angle[i]);
    eside[i] = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
    printf("%s;", str1);

    /* ge time spectra */

    sprintf(str1, "tge%3.3i", i);
    sprintf(str2, "ge det %3.3i time @ angle %6.2f", i, (float) angle[i]);
    tge[i] = mkTH1D(str1, str2, L14BITS, 1, M14BITS);
    printf("%s;", str1);

    /* ge low res energy spectra */

    sprintf(str1, "elo%3.3i", i);
    sprintf(str2, "det %3.3i ge low res  @ angle %6.2f", i, (float) angle[i]);
    elo[i] = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
    printf("%s;", str1);

    printf("\n");

    if (BinBGO)
    {

      /* ebgo spectra */

      sprintf(str1, "ebgo%3.3i", i);
      sprintf(str2, "det %3.3i BGO energy @ angle %6.2f", i, (float) angle[i]);
      ebgo[i] = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
      printf("%s;", str1);

      /* tbgo spectra */

      sprintf(str1, "tbgo%3.3i", i);
      sprintf(str2, "det %3.3i BGO time @ angle %6.2f", i, (float) angle[i]);
      tbgo[i] = mkTH1D(str1, str2, L14BITS, 1, M14BITS);
      printf("%s;", str1);

    };

    if (GeETmats)
    {

      /* create germanium et matrices */

      sprintf(str1, "geet_%3.3i", i);
      sprintf(str2, "det %3.3i ET @ angle %6.2f", i, (float) angle[i]);
      GeTEmat[i] = mkTH2F(str1, str2,
			  GeETmats_e_nn, GeETmats_e_lo, GeETmats_e_hi,
			  GeETmats_t_nn, GeETmats_t_lo, GeETmats_t_hi);
      printf("%s;", str1);
      GeTEmat[i]->SetXTitle("GS ge Energy");
      GeTEmat[i]->SetYTitle("GS ge Time");

    };


  };

  /* tac1 spectrum */

  sprintf(str1, "tac1");
  sprintf(str2, "tac1");
  tac1 = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
  printf("%s;", str1);

  /* sumehi spectrum */

  sprintf(str1, "sumehi");
  sprintf(str2, "sumehi");
  sumehi = mkTH1D(str1, str2, L14BITS, 1, M14BITS);
  printf("%s;", str1);

  /* sumelo spectrum */

  sprintf(str1, "sumelo");
  sprintf(str2, "sumelo");
  sumelo = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
  printf("%s;", str1);

  if (BinBGO)
  {

    /* sumeBGO spectrum */

    sprintf(str1, "sumebgo");
    sprintf(str2, "sumebgo");
    sumeBGO = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
    printf("%s;", str1);

    /* sumtbgo spectrum */

    sprintf(str1, "sumtbgo");
    sprintf(str2, "sumtbgo");
    sumtbgo = mkTH1D(str1, str2, L14BITS, 1, M14BITS);
    printf("%s;", str1);

    /* hitpattern */

    printf("\n");
    sprintf(str1, "bgohitpat");
    sprintf(str2, "Gammasphere bgo hitpattern");
    bgohitpat = mkTH1D(str1, str2, NGE, 0, NGE - 1);
    printf("%s;", str1);

    /* segment hit pattern */

    for (i = 0; i < NGE; i++)
    {
      gehits[i] = 0;
      for (j = 0; j < 8; j++)
	bgoseghit[i][j] = 0;
    };

  };

  /* sge spectrum (EFF sum) */

  sprintf(str1, "sge");
  sprintf(str2, "sge");
  sge = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
  printf("%s;", str1);

  /* sbgo spectrum (EFF sum) */

  sprintf(str1, "sbgo");
  sprintf(str2, "sbgo");
  sbgo = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
  printf("%s;", str1);

  /* sumeside spectrum */

  sprintf(str1, "sumeside");
  sprintf(str2, "sumeside");
  sumeside = mkTH1D(str1, str2, L12BITS, 1, M12BITS);
  printf("%s;", str1);

  /* sumtge spectrum */

  sprintf(str1, "sumtge");
  sprintf(str2, "sumtge");
  sumtge = mkTH1D(str1, str2, L14BITS, 1, M14BITS);
  printf("%s;", str1);

  /* hit patterns */

  printf("\n");
  sprintf(str1, "gehitpat");
  sprintf(str2, "Gammasphere ge hitpattern");
  gehitpat = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  printf("\n");
  sprintf(str1, "puhitpat");
  sprintf(str2, "Gammasphere pileup hitpattern");
  puhitpat = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  printf("\n");
  sprintf(str1, "orhitpat");
  sprintf(str2, "Gammasphere over-range hitpattern");
  orhitpat = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  printf("\n");
  sprintf(str1, "sidehitpat");
  sprintf(str2, "Gammasphere side chan ge hitpattern");
  sidehitpat = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  /* event length */

  sprintf(str1, "len");
  sprintf(str2, "event length");
  len = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  /* event length */

  sprintf(str1, "len_clean");
  sprintf(str2, "clean ge multiplicity");
  len_clean = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  /* clean+dirty multiplicity */

  sprintf(str1, "len_dirty");
  sprintf(str2, "clean+dirty ge multiplicity");
  len_dirty = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  /* BGO mult */

  sprintf(str1, "len_bgo");
  sprintf(str2, "BGO multiplicity");
  len_bgo = mkTH1D(str1, str2, NGE, 0, NGE - 1);
  printf("%s;", str1);

  /* ring spectra */

  printf("\n");
  for (i = 1; i <= 17; i++)
  {
    sprintf(str1, "r%3.3i", i);
    sprintf(str2, "ring %3.3i  @ angle %6.2f", i, (float) ring_angle[i]);
    ring[i] = mkTH1D(str1, str2, L14BITS, 1, M14BITS);
    printf("%s; ", str1);
  };
  printf("\n");

  /* other GS 2D spectra */

  printf("\n");
  sprintf(str1, "tmul_ehi");
  sprintf(str2, "tmul_ehi");
  TMul_ehi = mkTH2F(str1, str2, NGE, 0, NGE - 1, L11BITS, 1, M11BITS);
  printf("tmul_ehi; ");

  sprintf(str1, "dtbtev");
  sprintf(str2, "usec between events");
  dtbtev = mkTH1D(str1, str2, L10BITS, 1, M10BITS);

  sprintf(str1, "dtbtgm");
  sprintf(str2, "dt bt gammas");
  dtbtgm = mkTH1D(str1, str2, L14BITS, 1, M14BITS);
  printf("\n");

  /* make all gamma gamma matrices */

  for (i = 0; i < nGamxGam; i++)
  {
    printf("\n");
    ehi_ehi[i] = mkTH2F(GG[i].name, GG[i].name,
			GG[i].len, GG[i].lo, GG[i].hi,
			GG[i].len, GG[i].lo, GG[i].hi);
    printf("gamma gamma matrix #%i, ,<%s>\n", i, GG[i].name);
    printf("__%i ch (", GG[i].len);
    printf("%f -> %f)\n", GG[i].lo, GG[i].hi);

    if (GG[i].prompt_cond != -1)
      printf("(prompt) RF time mask #%i, <%s> [%s]\n",
	     GG[i].prompt_cond, TMask[GG[i].prompt_cond].name, TMask[GG[i].prompt_cond].strrange);
    else
      printf("no (prompt) RF time mask\n");

    if (GG[i].dtbggm_cond != -1)
      printf("delta time mask #%i, <%s> [%s]\n",
	     GG[i].dtbggm_cond, TMask[GG[i].dtbggm_cond].name, TMask[GG[i].dtbggm_cond].strrange);
    else
      printf("no delta time mask\n");

  }

  /* make the fera 1D spectra */

  if (NPEv1DBIN > 0)
  {
    printf("\n");
    printf("simple 1D pev binned spectra\n");
    printf("\n");

    for (i = 0; i < NPEv1DBIN; i++)
    {
      PEv1D[i] = mkTH1D(PEv[PB1D[i].pev].name, PEv[PB1D[i].pev].name,
		     PB1D[i].len, (Axis_t) PB1D[i].lo, (Axis_t) PB1D[i].hi);
      sprintf(str, "%s,pev[%i]", PEv[PB1D[i].pev].name, PB1D[i].pev);
      PEv1D[i]->SetXTitle(str);
      PEv1D[i]->SetYTitle("counts");
      printf("%s %ich; ", PEv[PB1D[i].pev].name, PB1D[i].len);
      printf("(%f->%f) [%p]\n", (float) PB1D[i].lo, (float) PB1D[i].hi, PEv1D[i]);
      fflush(stdout);
    };

  };


  /* make rate spectra */

  printf("nratesp=%i\n", nratesp);
  if (nratesp > 1)
  {

    printf("\n");
    printf("make rate spectra\n");
    printf("\n");

    for (i = 1; i < nratesp; i++)
    {
      sprintf(str, "%s_rate", ratesp[i].name);
      rate_sp[i] = mkTH1D(str, str,
	       ratesp[i].len, (Axis_t) ratesp[i].lo, (Axis_t) ratesp[i].hi);
      printf("made rate spectrum <%s>", str);
      sprintf(str, "%s,pev[%i]", ratesp[i].name, ratesp[i].pev);
      rate_sp[i]->SetXTitle(str);
      rate_sp[i]->SetYTitle("rate [Hz]");
      printf("__ncan: %i,", ratesp[i].len);
      printf("from %f ", (float) ratesp[i].lo);
      printf("to %f;", (float) ratesp[i].hi);
      printf("incrmt: %f", (float) ratesp[i].incrmt);
      printf("\n");

      fflush(stdout);
    };

  };

  /* make the fera 2D spectra */

  if (NPEv2DBIN > 0)
  {
    printf("\n");
    printf("simple 2D pevXpev binned spectra\n");
    printf("\n");
    for (i = 0; i < NPEv2DBIN; i++)
    {
      sprintf(str, "%sX%s", PEv[PB2D[i].pevx].name, PEv[PB2D[i].pevy].name);
      printf("%3i: %s\n", i, str);
      PEv2D[i] = mkTH2F(str, str,
			PB2D[i].lenx, PB2D[i].lox, PB2D[i].hix,
			PB2D[i].leny, PB2D[i].loy, PB2D[i].hiy);
      sprintf(str, "%s,pev[%i]", PEv[PB2D[i].pevx].name, PB2D[i].pevx);
      PEv2D[i]->SetXTitle(str);
      sprintf(str, "%s,pev[%i]", PEv[PB2D[i].pevy].name, PB2D[i].pevy);
      PEv2D[i]->SetYTitle(str);
      printf("__x: %i, %f-%f\n", PB2D[i].lenx, PB2D[i].lox, PB2D[i].hix);
      printf("__y: %i, %f-%f\n", PB2D[i].leny, PB2D[i].loy, PB2D[i].hiy);
      fflush(stdout);
    };

  };

  /* make the gamxpev 2D spectra */

  if (NGamPEv2DBin > 0)
  {
    printf("\n");
    printf("simple 2D gam x pev binned spectra\n");
    printf("\n");
    for (i = 0; i < NGamPEv2DBin; i++)
    {
      GamPEv2D[i] = mkTH2F(PEvG[i].name, PEvG[i].name,
			   PEvG[i].pevlen, PEvG[i].pevlo, PEvG[i].pevhi,
			   PEvG[i].gamlen, PEvG[i].gamlo, PEvG[i].gamhi);
      printf("%s \n", str);
      printf("__x: %i, %f-%f\n", PEvG[i].pevlen, PEvG[i].pevlo, PEvG[i].pevhi);
      printf("__y: %i, %f-%f\n", PEvG[i].gamlen, PEvG[i].gamlo, PEvG[i].gamhi);
      PEvG[i].prompt_cond = -1;
      fflush(stdout);
    };
  };

  printf("\n");
  printf("rate spectra\n");
  printf("\n");

  sprintf(str1, "evrate");
  sprintf(str2, "evrate, events/sec");
  evrate = mkTH1D(str1, str2, LENEVTIME, 1, LENEVTIME);
  sprintf(str1, "Seconds");
  evrate->SetXTitle(str1);
  evrate->SetYTitle("Hz");
  printf("evrate; ");

  sprintf(str1, "ehirate");
  sprintf(str2, "ehirate, events/sec");
  ehirate = mkTH1D(str1, str2, LENEVTIME, 1, LENEVTIME);
  sprintf(str1, "Seconds");
  ehirate->SetXTitle(str1);
  ehirate->SetYTitle("ehi Hz");
  printf("ehirate; ");

  sprintf(str1, "datarate");
  sprintf(str2, "datarate, bytes/sec");
  datarate = mkTH1D(str1, str2, LENDATATIME, 0, LENDATATIME - 1);
  sprintf(str1, "Seconds");
  datarate->SetXTitle(str1);
  datarate->SetYTitle("Hz");
  printf("datarate; ");
  printf("\n");
  printf("\n");

  /* make the special devmatrices */

  for (i = 0; i < nDevMatrix; i++)
  {
    i1 = DevMatrix[i].dhi - DevMatrix[i].dlo + 1;
    d1 = 0;
    d2 = (double) i1;
    sDevMatrix[i] = mkTH2F(DevMatrix[i].name, DevMatrix[i].name,
			   i1, d1, d2,
     DevMatrix[i].nn, (double) DevMatrix[i].ylo, (double) DevMatrix[i].yhi);
    printf("made devmatrix <%s>; ", DevMatrix[i].name);
    printf(" [%p]\n", sDevMatrix[i]);
    printf("__x %i ch from %f to %f; ", i1, (float) d1, (float) d2);
    printf("y %i ch from %f to %f\n", DevMatrix[i].nn, (float) DevMatrix[i].ylo, (float) DevMatrix[i].yhi);
  };
  printf("\n");

  /* make the 2D hit matrices */

  for (i = 0; i < nHit2DMat; i++)
  {
    i1 = Hit2DMat[i].dhi1 - Hit2DMat[i].dlo1 + 1;
    i2 = Hit2DMat[i].dhi2 - Hit2DMat[i].dlo2 + 1;
    d1 = 0;
    d2 = (double) i1 - 1;
    d3 = (double) i2 - 1;
    sHit2DMat[i] = mkTH2F(Hit2DMat[i].name, Hit2DMat[i].name,
			  i1, d1, d2,
			  i2, d1, d3);
    printf("made hit2dmat <%s>; ", Hit2DMat[i].name);
    printf(" [%p]\n", Hit2DMat[i]);
    printf("__x %i ch from %f to %f; ", i1, (float) d1, (float) d2);
    printf("y %i ch from %f to %f\n", i2, (float) d1, (float) d3);

  };
  printf("\n");

  /*------------------------*/
  /* execute user init code */
  /*------------------------*/

  printf("\n");
  printf("executing UserInit.h code\n");
  printf("\n");
#include "UserInit.h"

  /* update shared mem with minimal info       */
  /* so it is not empty before the first update */
  /* shared memory wellness checkpoint          */

  if (UseShareMemFile)
  {
    printf("\n");
    printf("Note: if the command below fails,\n");
    printf("increase the shared memory size!\n");
    printf("\n");
    printf("updating empty shared mem file... ");
    fflush(stdout);
    mfile->Update();
    printf("Done!\n");
    printf("\n");
    fflush(stdout);
  };

  /*----------------------*/
  /* setup signal catcher */
  /*----------------------*/

#ifdef LINUX
  signal(SIGHUP, signal_catcher);
#endif
#ifdef SOLARIS
  sigset(SIGHUP, signal_catcher);
#endif














  /*-------------------------------------------*/
  /* read data from tape and increment spectra */
  /*-------------------------------------------*/

  printf("\n");
  printf("sorting...\n");
  fflush(stdout);
  printf("\n");

  st = 0;
  CurFile = 0;
  CurEvNo = 0;
  NprintEvNo = 0;
  eov = 0;
  oldst = 0;
  tstart = time(NULL);
  tStart = time(NULL);
  tdmplast = time(NULL);
  while (st >= 0 && CurEvNo < nEvents && CurFile < Nfiles && eov == 0)
  {

    /*-----------------------------------*/
    /* mark pseudo event vector unfilled */
    /* and zero values (for modpev)..... */
    /*-----------------------------------*/

    for (k = 0; k < LenPEv; k++)
    {
      PEv[k].filled = FALSE;
      PEv[k].val = 0;
    };
    dclock1 = 0;
    dclock2 = 0;

    /*----------------*/
    /* get next event */
    /*----------------*/

    if (InputSrc == TAPE)
      st = get_gsII_ev(&hdr, ev, ext);	/* tape */
    else if (InputSrc == DISK)
      st = get_disk_ev(&hdr, ev, ext);	/* disk */
    else if (InputSrc == NET)
      st = get_net_ev(&hdr, ev, ext);	/* net */
      
#if(0)
    /*debug print clean bgo times */
    
    if (hdr.len_bgo>0 )
      {
#include "GSPrintEvent.h"
      }; 
#endif      

    /* update error/info statisics */

    if (st >= 0 && st < NERR)
      GSSortError[st]++;

    CurEvNo++;
    
    /* only process if we are past firstevent */
    
    if (CurEvNo < firstEvent)
      goto skipevent;
      
    NprintEvNo++;
    
    if (NprintEvNo <= NumToPrint)
    {
      printf("\n++++ CurEvNo: %3i ++++; ", CurEvNo);
      printf("buffer # %10i\n", GSSortError[1] + GSSortError[2] + GSSortError[3]);
    };

    /* if first event, store microscond clock */
    
    if(firstevent)
      if (st==0)
        {
        firstevent=0;
        start_ttL=hdr.ttL;
        start_ttM=hdr.ttM;
        start_ttH=hdr.ttH;
        };


    /* simple statics */

    NoCleanGE += hdr.len_clean;
    NoDirtyGE += hdr.len_dirty;
    NoCleanBGO += hdr.len_bgo;

    if (st == SUCCESS && hdr.len_clean > 0)
    {

      /* undo honeycomb suppression here */

      if (ignoreHS)
	for (j = 0; j < hdr.len_clean; j++)
	  if (ev[j].hs)
	  {
	    ev[j].bgohit = 0;
	    ev[j].hs = 0;
	  };

      /* honeycomb suppression statistics */

      for (j = 0; j < hdr.len_clean; j++)
	if (ev[j].hs)
	  NoCleanWithHS++;
	else if (ev[j].hs == 0 && ev[j].bgohit != 0)
	{
	  NoCleanWithHSmarked++;
	  NoCleanWithoutHS++;
	} else
	  NoCleanWithoutHS++;

    };




    /*----------------------------------------*/
    /* allow user to manipulate raw data here */
    /*----------------------------------------*/

#include "UserRawEv.h"




    /* houskeeping to do when we hit a buffer */

    if (st == GSIIERRNEWBUF)
    {

      /* data rate based on buffer reads */

      tnow = time(NULL);
      tnow -= tstart;
      d1 = tnow / 60.0;
      if (EvTime > 0 && EvTime < LENEVTIME)
	datarate->Fill(EvTime, DATASIZE);

      /* calc time since last dump */

      tdmp = time(NULL);
      tdmp -= tdmplast;
      tdmp /= 60;		/* now minutes */

      /* redefine event to be success if there is data */
      /* End of buffer reached does NOT necessarely    */
      /* mean the returned data is bad                 */

      if (ext[0] > 0 || hdr.len > 0)
      {
	st = SUCCESS;
	GSSortError[st]++;
      };

    };

    if (st < 0)
      printf("receiced st=%i at event no: %i\n", st, CurEvNo);

    /* print GS data */

    if (NprintEvNo <= NumToPrint)
      if (st > 0)
	printf("no data returned, return status: %i\n", st);


    /* interpret fera data */

    if (st == GOODEV && ext[0] > 0 && !DoNotInterpExt)
    {

      /* trap potentially bad external data */

      if (ext[0] >= (EXTLEN - 30))
	ext[0] = 0;

      /* print raw external detectors */

      if (NprintEvNo <= NumToPrint)
      {
	if (WriteRawFera)
	  if (ext[0] > 0)
	  {
	    printf("external data (%i raw words):\n", ext[0]);
	    fflush(stdout);
	    for (j = 1; j <= ext[0]; j++)
	    {

	      printf("%2i --> %10i, 0x%4.4x; ", j, ext[j], ext[j]);
	      fflush(stdout);

	      l = 32768;
	      printf("|");
	      for (k = 0; k < 16; k++)
	      {
		if ((ext[j] & l) == l)
		  printf("1");
		else
		  printf("0");
		l = l / 2;
		if ((k + 1) % 4 == 0)
		  printf("|");
	      };
	      printf("; ");
	      fflush(stdout);

	      printf("(");
	      printf("%1i,", (0x8000 & ext[j]) >> 15);
	      printf("%2i,", (0x7800 & ext[j]) >> 11);
	      printf("%4i", (0x07ff & ext[j]));
	      printf("),");

	      printf("\n");
	      fflush(stdout);
	    };

	  };
      };

      /* interpret external data */

      if (IsomerData)
	fera_st = decode_isomer(ext, &nEXT, fera_type,
				fera_vsn, fera_ch, fera_data, hdr.tac2);
      else
      fera_st = decode_fera(ext, &nEXT, fera_type,
			    fera_vsn, fera_ch, fera_data,
			    &dclock1, &dclock2);

      if (NprintEvNo <= NumToPrint)
	if (fera_st != SUCCESS)
	  printf("failed to interpret FERA data\n");

      /* sanitize ext data wrt binning */

      if (nEXT >= EXTLEN || fera_st != SUCCESS)
      {
	nEXT = 0;
	fera_st = -1;
      };

      for (i = 0; i < nEXT; i++)
      {
	if (fera_type[i] >= M12BITS)
	  fera_type[i] = M12BITS;
	if (fera_vsn[i] >= MAXVSN)
	  fera_vsn[i] = MAXVSN - 1;
	if (fera_ch[i] >= MAXFERACH)
	  fera_ch[i] = MAXFERACH - 1;
#if(0)
	if (fera_data[i] >= M12BITS)
	  fera_data[i] = M12BITS;
#endif
      };

    } else
    {
      nEXT = 0;
      fera_st = -1;
    };

    /* count files and check for eom */

    if (st == FILEMARK)
    {
      CurFile++;
      if (oldst == FILEMARK)
      {
	eov = 1;
	printf("\n");
	printf("found end-of-medium (double eof)\n");
	printf("\n");
      };
    };
    oldst = st;


    /*-----------------------------------------------------------*/
    /* dump all spectra on signal or dump every DumpEvery events */
    /*-----------------------------------------------------------*/

    if (WeWereSignalled || (int) tdmp >= DumpEvery)
    {

      /* disarm signal */

      WeWereSignalled = FALSE;

      /* check for command file */

      fp = fopen(CommandFileName, "r");
      if (fp != NULL)
      {

	printf("found command file: %s\n", CommandFileName);
	fgets(str, STRLEN, fp);
	printf("with command: %s\n", str);

	if ((p = strstr(str, "dumpevery")) != NULL)
	{

	  sscanf(str, "%s %i", str1, &DumpEvery);
	  printf("will dump to output file every %i minutes\n", DumpEvery);
	  fflush(stdout);

	} else if ((p = strstr(str, "printevents")) != NULL)
	{
	  /* reset print event counter */

	  nret = sscanf(str, "%s %i", str1, &i1);
	  if (nret == 2)
	    NumToPrint = i1;
	  NprintEvNo = 0;
	  SetFERAErrorPrint(NFERAErrToPrint);

	} else if ((p = strstr(str, "status")) != NULL)
	{
	  /* reset time 0 subtraction */

	  printf("current event time   is %f seconds\n", (float) EvTime);
	  printf("current event time 0 is %f seconds\n", (float) EvTime0);
	  r1 = EvTime + (float) EvTime0;	/* restore to true event time */
	  printf("true current event time is %f seconds\n", r1);
	  printf("Beta= %f\n", (float) Beta);
	  printf("CurEvNo= %i; ", CurEvNo);
	  printf("CurFile= %i\n",CurFile );
	  fflush(stdout);

	} else if ((p = strstr(str, "segmentHitPattern")) != NULL)
	{

	  printBGOSegmentHits();
	  fflush(stdout);

	} else if ((p = strstr(str, "resett0")) != NULL)
	{
	  /* reset time 0 subtraction */

	  printf("current event time 0 was %f seconds\n", (float) EvTime0);
	  printf("current event time is %f seconds\n", (float) EvTime);
	  r1 = EvTime + (float) EvTime0;	/* restore to true event time */
	  printf("true current event time is %f seconds\n", r1);
	  EvTime0 = r1;
	  printf("resetting t0 to %f seconds, true current event time\n", (float) EvTime0);

	  fflush(stdout);

	} else if ((p = strstr(str, "resetcounters")) != NULL)
	{
	  /* reset print event counter */

	  NprintEvNo = 0;
	  SetFERAErrorPrint(NFERAErrToPrint);

	  ResetGSII_bufpCounters();

	} else if ((p = strstr(str, "evtime0")) != NULL)
	{
	  nret = sscanf(str, "%s %f", str1, &r1);
	  CheckNoArgs(nret, 2, str);
	  EvTime0 = (double) r1;
	  printf("will start event time from %f seconds\n", (float) EvTime0);
	  fflush(stdout);

	} else if ((p = strstr(str, "stopsort")) != NULL)
	{
	  /* simulate end of tape to stop sort */

	  eov = 1;

	} else if ((p = strstr(str, "zapcounters")) != NULL)
	{

	  /* zap various counters so we */
	  /* can obtain a current estimate */
	  /* or various errors and such */

	  tstart = time(NULL);
	  tStart = time(NULL);
	  CurEvNo = 0;
	  firstEvent = 0;
	  nMultOK = 0;
	  NoCleanGE = 0;
	  NoDirtyGE = 0;
	  NoCleanBGO = 0;

	  for (i = 0; i < NERR; i++)
	    GSSortError[i] = 0;

	  ZapFERAErrors();

         /* segment hit patterns */

         for (i = 0; i < NGE; i++)
           {
           gehits[i] = 0;
           for (j = 0; j < 8; j++)
	   bgoseghit[i][j] = 0;
           };

	} else if ((p = strstr(str, "zapall")) != NULL)
	{

	  /* zap spectra */

	  zlist = gDirectory->GetList();
	  hiterator = zlist->MakeIterator();
	  while (htemp = (TH1 *) hiterator->Next())
	    htemp->Reset();
	  printf("all spectra were zapped @ ");
	  time_stamp();
	  fflush(stdout);

	  /* update */

	  if (!UseShareMemFile)
	  {
	    /* do nothing */
	  } else
	  {
	    printf("updating shared mem... ");
	    UPDSSHMEM
	  };

	} else if ((p = strstr(str, "beta")) != NULL)
	{
	  /* get beta value */

	  sscanf(str, "%s %f", str1, &Beta);
	  printf("setting beta to %9.3f\n", Beta);
	  SetBeta();

	} else if ((p = strstr(str, "zap")) != NULL)
	{
	  /* extract spectrum name */

	  sscanf(str, "%s %s", str1, spname);
	  htemp = (TH1D *) gROOT->FindObject(spname);

	  /* zap spectrup */

	  htemp->Reset();

	  printf("spectrum %s zapped @ ", spname);
	  time_stamp();
	  fflush(stdout);

	  /* update */

	  if (!UseShareMemFile)
	  {
	    /* do nothing */
	  } else
	  {
	    printf("updating shared mem... ");
	    UPDSSHMEM
	  };

	} else
	  printf("command \"%s\" not understood\n", str);

	/* delete command file */

	fclose(fp);
	sprintf(str, "\\rm %s", CommandFileName);
	system(str);
	printf("%s\n", str);

      } else
      {

	/* update sh mem or writeout root file */

	printf("time since last dump: %i minute(s)\n", (int) tdmp);
	tdmp = 0;
#if(0)
	printf("no command file: %s was found\n", CommandFileName);
#endif
        d1=start_ttL + 65536.0 * start_ttM + 2147483648.0 * start_ttH;
        d1/=1000000.0;
        runTime=(EvTime+EvTime0)-d1;
        
	WriteStatistics(CurEvNo, nMultOK, GSSortError,runTime);
	if (!UseShareMemFile)
	{

	  printf("*---------------------------------\n");
	  printf("* you cannot update a disk file.  \n");
	  printf("  you must wait for sort to finish\n");
	  printf("  or stop the sort! Ignoring you...\n");
	  printf("*---------------------------------\n");

	} else
	{
	  printf("updating shared mem... ");
	  UPDSSHMEM
	};

	tdmplast = time(NULL);

      };

      printf("continuing the sort...\n");
      fflush(stdout);

    };

    /* allow user to get first crack */
    /* at processing good event if they want to */
    /* now that FERAs have been interpreted. */
    /* The user can redefine st to st=GSIIREJECTEV */
    /* here and reject the event. That is how it */
    /* is done for isomer processing */

    if (st == GOODEV)
    {
#include "UserGoodEv.h"
    };

    /*---------------------*/
    /* process good events */
    /*---------------------*/

    if (st == GOODEV)
    {


      /* event time  */
      /* take v.r.t user specified t0 */

      EvTime = hdr.ttL + 65536.0 * hdr.ttM + 2147483648.0 * hdr.ttH;
      EvTime /= 1000000.0;	/* now seconds */
      EvTime -= EvTime0;

      /* update event rate spectrum */

      if (EvTime > 0 && EvTime < LENEVTIME)
	evrate->Fill(EvTime);

      /* time between events */

      d1 = (double) hdr.ttL - hdr_ttL_last;
      if (d1 < M10BITS && d1 >= 0)
	dtbtev->Fill(d1);
      hdr_ttL_last = hdr.ttL;

      /*-----------------------------------*/
      /* modify event with low energy cuts */
      /* count for multiplicity cuts       */
      /*-----------------------------------*/

      tmp_lc = hdr.len_clean;
      tmp_ld = hdr.len_dirty;
      tmp_lb = hdr.len_bgo;

      for (i = 0; i < hdr.len; i++)
      {

	if (ev[i].ehi > 0 && ev[i].ebgo == 0)
	{

	  /* clean ge */

	  if (ev[i].ehi <= egemin)
	  {
	    ev[i].id = 0;
	    ev[i].ehi = 0;
	    ev[i].elo = 0;
	    ev[i].eside = 0;
	    ev[i].tge = 0;
	    tmp_lc--;
	  };

	} else if (ev[i].ehi > 0 && ev[i].ebgo > 0)
	{

	  /* dirty ge */

	  if (ev[i].ehi < egemin)
	  {
	    ev[i].id = 0;
	    ev[i].ehi = 0;
	    ev[i].elo = 0;
	    ev[i].eside = 0;
	    ev[i].tge = 0;
	  };

	  if (ev[i].ebgo < ebgomin)
	  {
	    ev[i].ebgo = 0;
	    ev[i].tbgo = 0;
	  };

	  /* if both are zero after above eval, */
	  /* then we lost a dirty germanium     */

	  if (ev[i].ehi == 0 && ev[i].ebgo == 0)
	    tmp_ld--;

	} else if (ev[i].ehi == 0 && ev[i].ebgo > 0)
	{

	  /* clean bgo */

	  if (ev[i].ebgo < ebgomin)
	  {
	    ev[i].ebgo = 0;
	    ev[i].tbgo = 0;
	    tmp_lb--;
	  };

	};

      };

      /* find 'real' multiplicities for event processing */

      lenclean = tmp_lc;
      lendirty = lenclean + tmp_ld;
      lentotal = lendirty + tmp_lb;

      /* modify multiplicities for hs */

      for (j = 0; j < hdr.len_clean; j++)
	if (ev[j].hs)
	{
	  lenclean--;		/* true lenclean */
	  lendirty++;		/* true lendirty */
	};

      /* check whether event fulfills multiplicity requirements */
      /* (affects ALL binning!)                                 */

      MultOK = 0;
      if (lenclean >= CleanMultLo && lenclean <= CleanMultHi)
	if (lendirty >= DirtyMultLo && lendirty <= DirtyMultHi)
	  if (lentotal >= TotalMultLo && lentotal <= TotalMultHi)
	    MultOK = 1;

      /* debug print multiplicity test results */

      if (NprintEvNo <= NumToPrint)
      {

	printf("event time %20.6f (with t0=%20.6f) seconds; ~", (float) EvTime, (float) EvTime0);

	i0 = (int) EvTime;
	i1 = i0 / 86400;
	i0 = i0 - i1 * 86400;
	i2 = i0 / 3600;
	i0 = i0 - i2 * 3600;
	i3 = i0 / 60;
	i0 = i0 - i3 * 60;
	printf("%3id", i1);
	printf("%2ih", i2);
	printf("%2im", i3);
	printf("%2is\n", i0);

	printf("clean: %i, ", lenclean);
	printf("clean+dirty: %i, ", lendirty);
	printf("total: %i; ", lentotal);
	fflush(stdout);

	if (MultOK)
	  printf("multiplicity requirements ok\n");
	else
	  printf("failed multiplicity requirements\n");

	fflush(stdout);

      };

      /* start binning if multiplicities were ok */

      if (MultOK)
      {

	/* count these events */

	nMultOK++;

	/* doppler correct the energies */
	/* and limit range              */

	for (j = 0; j < hdr.len; j++)
	{

	  /* hi res */

	  r1 = (float) ev[j].ehi + (float) rand() / RAND_MAX - 0.5;
	  if (ehiDoGainCor)
	    r1 = r1 * ehiGeGain[ev[j].id] + ehiGeOffset[ev[j].id];
	  r1 *= DopCorFacHires[ev[j].id];
	  ev[j].ehi = (unsigned short int) (r1 + 0.5);
	  if (ev[j].ehi >= M14BITS)
	    ev[j].ehi = 0;

	  /* lores */

	  r1 = (float) ev[j].elo + (float) rand() / RAND_MAX - 0.5;
	  if (eloDoGainCor)
	    r1 = r1 * eloGeGain[ev[j].id] + eloGeOffset[ev[j].id];
	  r1 *= DopCorFac[ev[j].id];
	  ev[j].elo = (unsigned short int) (r1 + 0.5);
	  if (ev[j].elo >= M14BITS)
	    ev[j].elo = 0;

	  /* eside */

	  r1 = (float) ev[j].eside + (float) rand() / RAND_MAX - 0.5;
	  if (esideDoGainCor)
	    r1 = r1 * esideGeGain[ev[j].id] + esideGeOffset[ev[j].id];
	  r1 *= DopCorFac[ev[j].id];
	  ev[j].eside = (unsigned short int) (r1 + 0.5);
	  if (ev[j].eside >= M14BITS)
	    ev[j].eside = 0;

	  /* keep an eye on ehi rate by itself */

	  if (EvTime > 0 && EvTime < LENEVTIME)
	    if (ev[j].ehi >= egemin)
	      ehirate->Fill(EvTime);

	  /* lo res */

	  r1 = (float) ev[j].elo + (float) rand() / RAND_MAX - 0.5;
	  r1 *= DopCorFac[ev[j].id];
	  ev[j].elo = (unsigned short int) (r1 + 0.5);
	  if (ev[j].elo >= M12BITS)
	    ev[j].elo = 0;

	  /* eside */

	  r1 = (float) ev[j].eside + (float) rand() / RAND_MAX - 0.5;
	  r1 *= DopCorFac[ev[j].id];
	  ev[j].eside = (unsigned short int) (r1 + 0.5);
	  if (ev[j].eside >= M12BITS)
	    ev[j].eside = 0;

	  if (BinBGO)
	  {

	    /* ebgo */

	    r1 = (float) ev[j].ebgo + (float) rand() / RAND_MAX - 0.5;
	    r1 *= DopCorFac[ev[j].id];
	    ev[j].ebgo = (unsigned short int) (r1 + 0.5);
	    if (ev[j].ebgo >= M12BITS)
	      ev[j].ebgo = 0;

	    /* limit data range for time spectra */

	    if (ev[j].tge >= M14BITS)
	      ev[j].tge = 0;
	    if (ev[j].tbgo >= M14BITS)
	      ev[j].tbgo = 0;

	  };

	};

	/* constrict other data ranges */
	/* this will be done in lower  */
	/* level routines later        */

	if (hdr.len >= NGE)
	  hdr.len = NGE - 1;
	if (hdr.len_clean >= NGE)
	  hdr.len_clean = 0;
	if (hdr.len_dirty >= NGE)
	  hdr.len_dirty = 0;
	if (hdr.len_bgo >= NGE)
	  hdr.len_bgo = 0;

	if (hdr.tac1 >= M12BITS)
	  hdr.tac1 = 0;
	if (hdr.tac2 >= M12BITS)
	  hdr.tac2 = 0;

	if (hdr.sumge >= M12BITS)
	  hdr.sumge = 0;
	if (hdr.sumbgo >= M12BITS)
	  hdr.sumbgo = 0;

	/*--------------------*/
	/* time manipulations */
	/*--------------------*/


	/* undo tac2 subtraction */
	/* (Note: BGO times are NOT subtracted in EFFs!) */

	if (undotac2sub)
	  for (j = 0; j < hdr.len; j++)
	    ev[j].tge += (int) ((float) hdr.tac2 * undotac2sub_factor
				- undotac2sub_offset);

	/* subtract tac2 rf from tge & tbgo, sanitize result */

	for (j = 0; j < hdr.len; j++)
	{

	  /* ge time */

	  if (ev[j].tge > 0)
	  {

	    tgerf[j] = (int) ev[j].tge - tGEoffset[ev[j].id];

	    if (subtac2GE)
	      tgerf[j] -= (int) ((float) hdr.tac2 * tac2FactorGE) + geRFOffset;

	    if (tgerf[j] < 0)
	      tgerf[j] = 0;
	    if (tgerf[j] >= M14BITS)
	      tgerf[j] = 0;
	  } else
	    tgerf[j] = 0;

	  /* bgo time */

	  if (ev[j].tbgo > 0)
	  {

	    tbgorf[j] = (int) ev[j].tbgo - tBGOoffset[ev[j].id];

	    /* manipulate? */

	    if (subtac2BGO)
	      tbgorf[j] -= (int) ((float) hdr.tac2 * tac2FactorBGO) + BGORFOffset;

	    /* sanitize */

	    if (tbgorf[j] < 0)
	      tbgorf[j] = 0;
	    if (tbgorf[j] >= M14BITS)
	      tbgorf[j] = 0;
	  } else
	    tbgorf[j] = 0;
	};

	/* print events */

	if (NprintEvNo <= NumToPrint)
	{
#include "GSPrintEvent.h"
	};

	/*--------------------------------------*/
	/* fill GammaSphere pseudo event vector */
	/*--------------------------------------*/

	PEv[GAMMASPHERE].val = 1;
	PEv[GAMMASPHERE].filled = TRUE;
	PEv[LEN_CLEAN].val = hdr.len_clean;
	PEv[LEN_CLEAN].filled = TRUE;
	PEv[LEN_DIRTY].val = hdr.len_dirty;
	PEv[LEN_DIRTY].filled = TRUE;
	PEv[LEN_BGO].val = hdr.len_bgo;
	PEv[LEN_BGO].filled = TRUE;

	PEv[LENCLEAN].val = lenclean;
	PEv[LENCLEAN].filled = TRUE;
	PEv[LENDIRTY].val = lendirty;
	PEv[LENDIRTY].filled = TRUE;
	PEv[TAC2].val = hdr.tac2;
	PEv[TAC2].filled = TRUE;
	PEv[LENEXT].val = ext[0];
	PEv[LENEXT].filled = TRUE;

	/*----*/
	/* HK */
	/*----*/

	if (mkHK == FALSE)
	{

	  /* just use what was calculated in the EFFs */

	  PEv[SUMGE].val = hdr.sumge;
	  PEv[SUMGE].filled = TRUE;
	  PEv[SUMBGO].val = hdr.sumbgo;
	  PEv[SUMBGO].filled = TRUE;
	  PEv[LENTOTAL].filled = TRUE;
	  PEv[LENTOTAL].val = lentotal;

	} else
	{

	  /* make our own calculation here in GSSort */

	  PEv[SUMGE].filled = TRUE;
	  PEv[SUMGE].val = 0;
	  PEv[SUMBGO].filled = TRUE;
	  PEv[SUMBGO].val = 0;
	  PEv[LENTOTAL].filled = TRUE;
	  PEv[LENTOTAL].val = 0;

	  i1 = 0;
	  i2 = hdr.len_clean;
	  for (i = i1; i < i2; i++)
	    if (tgerf[i] >= HKgelo && tgerf[i] <= HKgehi)
	    {
	      PEv[SUMGE].val += ev[i].elo;
	      PEv[LENTOTAL].val++;
	    };

	  i1 = hdr.len_clean;
	  i2 = i1 + hdr.len_dirty;
	  if (i2 > hdr.len)
	    i2 = hdr.len;
	  for (i = i1; i < i2; i++)
	  {
	    ok1 = (tgerf[i] >= HKgelo && tgerf[i] <= HKgehi);
	    ok2 = (tbgorf[i] >= HKBGOlo && tbgorf[i] <= HKBGOhi);
	    if (ok1)
	      PEv[SUMGE].val += ev[i].elo;
	    if (ok2)
	      PEv[SUMBGO].val += ev[i].ebgo;
	    if (ok1 || ok2)
	      PEv[LENTOTAL].val++;	/* only once */
	  };

	  i1 = hdr.len_clean + hdr.len_dirty;
	  i2 = i1 + hdr.len_bgo;
	  if (i2 > hdr.len)
	    i2 = hdr.len;
	  for (i = i1; i < i2; i++)
	    if (tbgorf[i] >= HKBGOlo && tbgorf[i] <= HKBGOhi)
	    {
	      PEv[SUMBGO].val += ev[i].ebgo;
	      PEv[LENTOTAL].val++;
	    };

	  hdr.sumbgo = (unsigned short int) (PEv[SUMBGO].val + 0.5);
	  hdr.sumge = (unsigned short int) (PEv[SUMGE].val + 0.5);


	};

	/* fill clock pevs */

	if (dclock1 > 0)
	{
	  PEv[CLOCK1].val = dclock1;
	  PEv[CLOCK1].filled = TRUE;
	};
	if (dclock2 > 0)
	{
	  PEv[CLOCK2].val = dclock2;
	  PEv[CLOCK2].filled = TRUE;
	};

	/* gamma gamma matrix update start status */

	for (i = 0; i < nGamxGam; i++)
	  GG[i].filled = TRUE;

	/* honeycomb suppessed clean germanium signals */

	j = 0;
	for (i = 0; i < hdr.len_clean && i < NGSPEV; i++)
	  if (!ev[i].hs)
	  {

	    PEv[j + GSPEVL1].val = (double) ev[i].ehi;
	    PEv[j + GSPEVL1].filled = TRUE;

	    PEv[j + GSPEVL2].val = (double) tgerf[i];
	    PEv[j + GSPEVL2].filled = TRUE;

	    j++;
	  };

	/*-------------------------------*/
	/* fill FERA pseudo event vector */
	/*-------------------------------*/

	/* __with threshold and energy calibration */

	for (j = 0; j < nEXT; j++)
	{
	  cpev = FERAPEvlookup[fera_vsn[j]][fera_ch[j]];
#if(0)
	  printf("\n");
	  printf("FERA wd # %2i; ", j);
	  printf("VSN=%3i; ", fera_vsn[j]);
	  printf("CH=%2i; ", fera_ch[j]);
	  printf("DATA=%5i; ", fera_data[j]);
	  printf("PEV #=%3i\n", cpev);

	  for (i = 0; i < MAXVSN; i++)
	    for (j = 0; j < MAXFERACH; j++)
	      if (FERAPEvlookup[i][j] > 0)
		printf("FERAPEvlookup[%i][%i]=%i\n", i, j, FERAPEvlookup[i][j]);
#endif
	  if (cpev > 0)
	    if (fera_data[j] >= PEv[cpev].thresh)
	    {
	      PEv[cpev].filled = TRUE;
	      if (PEv[cpev].calib)
	      {
		rn = (double) rand() / RAND_MAX - 0.5;
		PEv[cpev].val = (fera_data[j] + rn) * PEv[cpev].gain + PEv[cpev].offset;
	      } else
		PEv[cpev].val = fera_data[j];
	    };
	};

	/* -------------------------------------- */
	/* fill the special DSSD max energy pevs */
	/* (this is the same loop as I use later */
	/* to update the 2Dhit matrix, but un-  */
	/* fortunately I need to run it here    */
	/* first to find the energy and stuff   */
	/* it into a pev so it can be prosessed) */
	/* -------------------------------------- */

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

	  /* hunt for hits */

	  emax = 0;
	  for (j = Hit2DMat[i].dlo1; j <= Hit2DMat[i].dhi1; j++)
	    for (k = Hit2DMat[i].dlo2; k <= Hit2DMat[i].dhi2; k++)
	      if (PEv[j].filled && PEv[k].filled)
		if (PEv[j].val >= Hit2DMat[i].thresh)
		  if (PEv[k].val >= Hit2DMat[i].thresh)
		  {

		    /* and keep record of max energy */

		    if (PEv[j].val > emax)
		      emax = PEv[j].val;
		    if (PEv[k].val > emax)
		      emax = PEv[k].val;
#if(0)
		    printf("hit: %i, ", j - Hit2DMat[i].dlo1);
		    printf("%i, ", k - Hit2DMat[i].dlo2);
		    printf("%f, %f ", (float) PEv[j].val, (float) PEv[k].val);
		    printf("emax=%f\n", (float) emax);
#endif

		  };

	  /* update device (usually dsssd) energy */

	  if (emax > 0)
	  {
	    PEv[Hit2DMat[i].emaxpev].val = emax;
	    PEv[Hit2DMat[i].emaxpev].filled = TRUE;
#if(0)
	    printf("set pev# %i value: %f\n", Hit2DMat[i].emaxpev, PEv[Hit2DMat[i].emaxpev].val);
#endif
	  };

	};

	/* allow user access before gating */
	/* add your own pevs here */

#include "UserPreCond.h"


	/*------------------------------------------------*/
	/* do early/primary gating                        */
	/* typically energies gated with associated times */
	/*------------------------------------------------*/

	DoPrimaryGating();

	/*---------------------------------*/
	/* condition/counting/mpev ops etc */
	/*---------------------------------*/

	/* manipulate pseudo event vectors */

	mpevOP_reg();

	/* find conditions fulfilled */

	FindConditions();

	/* count conditions */

	CountConditions();

	/* find conditions fulfilled again, now we have counted */
	/* them there may be conditions on a counting pev       */

	FindConditions();

	/* evaluate the pesky counting pevs a few        */
	/* generations to make sure we catch ALL of them */

	for (k = 0; k < CountPevLevCheck; k++)
	{
	  mpevOP_special();
	  FindConditions();
	};

	/* print pev before conditions are applied */

	if (NprintEvNo <= NumToPrint && 0 && PrintRawPEv)
	{
	  printf("pevs before conditions:\n");
	  for (k = 0; k < LenPEv; k++)
	    if (PEv[k].filled)
	      printf("pev[%i]=%f,  %s\n", k, (float) PEv[k].val, PEv[k].name);
	};

	ApplyAllConditions();

	/* special treatment of GS pevs */

	if (!PEv[GAMMASPHERE].filled)
	{
	  PEv[LEN_CLEAN].filled = FALSE;
	  PEv[LEN_DIRTY].filled = FALSE;
	  PEv[LEN_BGO].filled = FALSE;
	  PEv[LENCLEAN].filled = FALSE;
	  PEv[LENDIRTY].filled = FALSE;
	  PEv[LENTOTAL].filled = FALSE;
	  PEv[SUMGE].filled = FALSE;
	  PEv[SUMBGO].filled = FALSE;
	  PEv[TAC2].filled = FALSE;
	  PEv[LENEXT].filled = FALSE;
	};

	/*---------------------------------------------*/
	/* all done filling/manipulating/counting pevs */
	/*---------------------------------------------*/

	/* print some events */

	if (NprintEvNo <= NumToPrint)
	{

	  /* print conditions */

	  printf("conditions ok: ");
	  for (k = 0; k < NPEvCond; k++)
	    if (PEvCond[k].ok)
	      printf("%s,", PEvCond[k].name);
	  printf("\n");

	  /* print pseudo event vector */

	  printf("final pev-->\n");
	  for (k = 0; k < LenPEv; k++)
	    if (PEv[k].filled)
	      printf("pev[%s,%i]=%f\n", PEv[k].name, k, PEv[k].val);
	  printf("\n");

	  /* print matrices */

	  printf("gamgam matrix status:\n");
	  for (k = 0; k < nGamxGam; k++)
	  {
	    printf("%2i: GG name [%s] ", k, GG[k].name);
	    printf("filled: %i\n", GG[k].filled);
	  };

	};

	/* increment clean spectra, time spectra, etc */

	if (PEv[GAMMASPHERE].filled)
	  for (j = 0; j < hdr.len_clean; j++)
	    if (ev[j].id >= 1)
	    {

	      /* __only update if not honeycomb suppressed */
	      /* and otherwise have a good event */

	      if (!ev[j].hs && !ev[j].over && !ev[j].pu)
	      {
		if (ev[j].ehi >= egemin)
		  ehi[ev[j].id]->Fill(ev[j].ehi, ACFac[ev[j].id]);
		elo[ev[j].id]->Fill(ev[j].elo, ACFac[ev[j].id]);
		if (ev[j].eside >= esidemin)
		  eside[ev[j].id]->Fill(ev[j].eside, ACFac[ev[j].id]);
		ring[angno[ev[j].id]]->Fill(ev[j].ehi, ACFac[ev[j].id]);

		sumehi->Fill(ev[j].ehi, ACFac[ev[j].id]);
		sumelo->Fill(ev[j].elo, ACFac[ev[j].id]);
		if (ev[j].eside >= esidemin)
		  sumeside->Fill(ev[j].eside, ACFac[ev[j].id]);

		if (GeETmats)
		  if (ev[j].ehi > GeETmats_e_lo && ev[j].ehi < GeETmats_e_hi)
		    if (ev[j].tge > GeETmats_t_lo && ev[j].tge < GeETmats_t_hi)
		      GeTEmat[ev[j].id]->Fill(ev[j].ehi, ev[j].tge);

		/* mult-gamma matrix */

		if (ev[j].ehi < M11BITS && lentotal < NGE)
		  TMul_ehi->Fill(lentotal, ev[j].ehi);

		/* update gamma-gamma matrices */
		/* __update controlled by GG[l].filled */

		for (l = 0; l < nGamxGam; l++)
		  if (GG[l].filled)
		    for (k = j + 1; k < hdr.len_clean; k++)
		      if (!ev[k].hs)
			if (ev[k].id >= 1)
			  if (ev[k].id <= MAXID)
			    if (ev[j].ehi < GG[l].hi && ev[j].ehi >= GG[l].lo)
			      if (ev[k].ehi < GG[l].hi && ev[k].ehi >= GG[l].lo)
			      {
				/* require valid promt time mask lookup */

				ok = TRUE;
				if (GG[l].prompt_cond != NOTDEF)
				{
				  ok = TMask[GG[l].prompt_cond].lot[tgerf[j]];
				  ok &= TMask[GG[l].prompt_cond].lot[tgerf[k]];
				};

				/* require valid dt mask lookup */

				if (ok && GG[l].dtbggm_cond != NOTDEF)
				{
				  i1 = tgerf[k] - tgerf[j];
				  if (i1 < 0)
				    i1 = -i1;
				  if (i1 >= L14BITS)
				    i1 = 0;
				  ok &= TMask[GG[l].dtbggm_cond].lot[i1];
				};

				/* update */

				if (ok)
				{
				  ehi_ehi[l]->Fill(ev[j].ehi, ev[k].ehi);
				  ehi_ehi[l]->Fill(ev[k].ehi, ev[j].ehi);
				};

			      };


	      };

	      /* update hitpatterns */

	      if (ev[j].ehi > egemin)
	      {
		gehitpat->Fill((double) ev[j].id, 1);
		gehits[ev[j].id]++;
	      };
	      if (ev[j].eside > esidemin)
		sidehitpat->Fill((double) ev[j].id, 1);
	      if (ev[j].pu)
		puhitpat->Fill((double) ev[j].id, 1);
	      if (ev[j].over)
		orhitpat->Fill((double) ev[j].id, 1);

	      /* notice, outside time cuts */

	      tge[ev[j].id]->Fill(tgerf[j], 1);
	      sumtge->Fill(tgerf[j], 1);

	    };

	/* time between gammas array */

	for (j = 0; j < hdr.len_clean; j++)
	  if (!ev[j].hs)
	    for (k = j + 1; k < hdr.len_clean; k++)
	      if (!ev[k].hs)
	      {
		i1 = tgerf[k] - tgerf[j];
		if (i1 < 0)
		  i1 = -i1;
		if (i1 >= 0 && i1 < M12BITS)
		  dtbtgm->Fill((float) i1, 1.0);
	      };
#if(0)
	nxx++;
	if (nxx > 20)
	  exit(0);
#endif

	if (PEv[GAMMASPHERE].filled)
	{


	  if (BinBGO)
	  {

	    /* bin BGOs */

	    for (j = hdr.len_clean; j < hdr.len; j++)
	      if (ev[j].id > 0 && ev[j].id <= MAXID)
	      {
		ebgo[ev[j].id]->Fill(ev[j].ebgo, ACFac[ev[j].id]);
		sumeBGO->Fill(ev[j].ebgo, ACFac[ev[j].id]);
		tbgo[ev[j].id]->Fill(tbgorf[j], 1);
		sumtbgo->Fill(tbgorf[j], 1);
		bgohitpat->Fill((double) ev[j].id, 1);

		/* segment hits (not a spectrum) */

		for (k = 0; k < 8; k++)
		  if (bincode[k] & ev[j].bgohit)
		    bgoseghit[ev[j].id][k]++;

	      };

	  };

	  /* other event information */

	  tac1->Fill(hdr.tac1, 1);
	  sge->Fill(hdr.sumge, 1);
	  sbgo->Fill(hdr.sumbgo, 1);

	  len->Fill(hdr.len, 1);
	  len_clean->Fill(hdr.len_clean, 1);
	  if ((hdr.len_clean + hdr.len_dirty) < NGE)
	    len_dirty->Fill(hdr.len_clean + hdr.len_dirty, 1);
	  len_bgo->Fill(hdr.len_bgo, 1);

	};

	/* increment pev 1D spectra */

	for (i = 0; i < NPEv1DBIN; i++)
	  if (PEv[PB1D[i].pev].filled)
	    if (PEv[PB1D[i].pev].val > PB1D[i].lo)
	      if (PEv[PB1D[i].pev].val < PB1D[i].hi)
		PEv1D[i]->Fill(PEv[PB1D[i].pev].val, 1);

	/* increment pev 2D spectra */

	for (i = 0; i < NPEv2DBIN; i++)
	  if (PEv[PB2D[i].pevx].filled)
	    if (PEv[PB2D[i].pevy].filled)
	      if (PEv[PB2D[i].pevx].val > PB2D[i].lox)
		if (PEv[PB2D[i].pevx].val < PB2D[i].hix)
		  if (PEv[PB2D[i].pevy].val > PB2D[i].loy)
		    if (PEv[PB2D[i].pevy].val < PB2D[i].hiy)
		      PEv2D[i]->Fill(PEv[PB2D[i].pevx].val, PEv[PB2D[i].pevy].val, 1);

	/* increment gam pev 2D spectra */

	if (PEv[GAMMASPHERE].filled)
	  for (i = 0; i < NGamPEv2DBin; i++)
	    if (PEv[PEvG[i].pev].filled)
	      if (PEv[PEvG[i].pev].val < PEvG[i].pevhi)
		if (PEv[PEvG[i].pev].val > PEvG[i].pevlo)
		  for (j = 0; j < hdr.len_clean; j++)
		    if (ev[j].ehi > PEvG[i].gamlo)
		      if (ev[j].ehi < PEvG[i].gamhi)
		      {
			ok = TRUE;
			if (PEvG[i].prompt_cond != -1)
			  ok &= TMask[PEvG[i].prompt_cond].lot[tgerf[j]];
			if (ok)
			  GamPEv2D[i]->Fill(PEv[PEvG[i].pev].val, ev[j].ehi, 1);
		      };

	/* increment devmatrices */

	for (i = 0; i < nDevMatrix; i++)
	  for (j = DevMatrix[i].dlo; j <= DevMatrix[i].dhi; j++)
	    if (PEv[j].filled)
	      if (PEv[j].val >= DevMatrix[i].ylo)
		if (PEv[j].val <= DevMatrix[i].yhi)
		  sDevMatrix[i]->Fill(j - DevMatrix[i].dlo, PEv[j].val, 1);

	/* increment hit2dmatrices  */

	for (i = 0; i < nHit2DMat; i++)
	  for (j = Hit2DMat[i].dlo1; j <= Hit2DMat[i].dhi1; j++)
	    for (k = Hit2DMat[i].dlo2; k <= Hit2DMat[i].dhi2; k++)
	      if (PEv[j].filled && PEv[k].filled)
		if (PEv[j].val >= Hit2DMat[i].thresh)
		  if (PEv[k].val >= Hit2DMat[i].thresh)
		    sHit2DMat[i]->Fill(j - Hit2DMat[i].dlo1, k - Hit2DMat[i].dlo2, 1);

	/* update any rate spectra */

	for (i = 0; i < LenPEv; i++)
	  if (PEv[i].ratespno > 0)
	    if (PEv[i].filled)
	      if (EvTime > ratesp[PEv[i].ratespno].lo && EvTime < ratesp[PEv[i].ratespno].hi)
		rate_sp[PEv[i].ratespno]->Fill(EvTime, ratesp[PEv[i].ratespno].incrmt);

	/*-------------------------*/
	/* execute user event code */
	/*-------------------------*/

#include "UserEv.h"

	/*------------------*/
	/* skip event point */
	/*------------------*/

    skipevent:
	{
	};

      };			/* if (MultOK) */

    };				/* if (st == GOODEV) */
    

  };				

  printf("\n");
  printf("Done sorting!\n");
  printf("\n");














  /*-------------*/
  /* close input */
  /*-------------*/

  if (InputSrc == NET)
    CloseUdpReceiver();

  /*-----------------------*/
  /* list the root spectra */
  /*-----------------------*/

#if(0)
  if (!UseShareMemFile)
  {
    printf("\n");
    printf("list of ROOT spectra:\n");
    printf("\n");
    TList          *hlist = gDirectory->GetList();
    hlist->ls();
    fflush(stdout);
  };
#endif

  /*------------*/
  /* statistics */
  /*------------*/

  d1=start_ttL + 65536.0 * start_ttM + 2147483648.0 * start_ttH;
  d1/=1000000.0;
  runTime=(EvTime+EvTime0)-d1;

  WriteStatistics(CurEvNo, nMultOK, GSSortError, runTime);

  /*-----------------------*/
  /* save all ROOT spectra */
  /*-----------------------*/

  if (!UseShareMemFile)
  {
    printf("attempting to close root file...");
    fflush(stdout);
    if (UpdateRootFile)
    {

      /* simple, we allready have file, */
      /* so just write to it */

      f1->Write();

    } else
    {

      /* go through all the crap */
      /* necessary to do this */

      WRITEALLHISTS

    }
    f1->Print();
    f1->Close();
    printf("done!\n");
    fflush(stdout);
  } else
  {
    UPDSSHMEM
      mfile->Print();
    printf("\n");
    mfile->ls();
    printf("\n");
  };
  printf("\n");

  /*----------------------*/
  /* condition statistics */
  /*----------------------*/

  printf("condition statistics:");
  printf("\n");

  for (i = 0; i < NPEvCond; i++)
    if (PEvCond[i].filled)
    {
      printf("condition #%3i: %9i, ", i, PEvCond[i].count);
      printf("%6.3f%%, name=%s\n", PEvCond[i].count * 100.0 / (float) GSSortError[0], PEvCond[i].name);
    };
  printf("\n");


  /* isomer statistics */

  if (IsomerData)
    wrIsomerStatistics(CurEvNo);

  /*-------------------------*/
  /* execute user exit code */
  /*-------------------------*/

  printf("\n");
  printf("executing UserExit.h code\n");
  printf("\n");
#include "UserExit.h"


  /*------*/
  /* done */
  /*------*/

  printf("\n");
  printf("$Id: GSSort.cxx,v 1.330 2011/03/07 20:09:10 tl Exp $\n");
  printf("*** all done on: ");
  time_stamp();
  printf("\n");
  return 0;

}

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

int
GSSort_read_chat(char *name)

{

  /* declarations */

  FILE           *fp, *fp1;
  char           *pc, *pc1, str[STRLEN] = {'0'}, str1[STRLEN] = {'0'}, str2[STRLEN] = {'0'};
  char            str3[STRLEN], str4[STRLEN], str5[STRLEN], str6[STRLEN];
  int             nn = 0, nni = 0, st, PType;
  char           *p;
  int             i, k, i1, i2, i3, i4, i5, i6;
  int             j1, j2, j3, j4, j5, j6, j7;
  float           f1, f2, f3, f4;
  int             echo = 0, nret;
  double          d1;

  /* prototypes */

  int             set_gsII_indrive(char *);
  TCutG          *rd2dwin(Char_t *);
  int             FindPEvMatNo(char *, int *);
  void            FindCondNo(char *, int *);
  int             SetFERAVSN(char *, int);
  void            InitializeFERAvsntable();
  void            ZeroFERAtypeVSN();
  void            PrintFERATypes();
  void            SetNPosWarn(int);
  void            SetRecordVer_tape(int);
  void            SetRecordVer_disk(int);
  void            SetRecordVer_net(int);
  int             str_decomp(char *, int, int *);
  void            FindTimeMaskNo(char *, int *);
  int             RdOffFile(char *, int *);
  int             RdGeCalFile(char *, float *, float *);
  void            CheckNoArgs(int, int, char *);
  void            SetSpecial(char *str);
  void            SetExportModNWords(char *, int);
  void            SetlongRangeTDCNWords(char *, int);
  void            setIsomerIDs(int);
  void            SetClockPar(int, float, float);
  void            SetFERAvsntable(int, int);
  void            SetnFeraDebug(int);

  /* initialize */

  if (!PevsInitialized)
  {
    /* make sure we only do this once          */
    /* since other chat files can be included  */
    /* and this function is called recursively */

    PevsInitialized = TRUE;

    for (i = 0; i < MAXPEV; i++)
      PEvCond[i].filled = FALSE;

    InitializeFERAvsntable();
    ZeroFERAtypeVSN();

    for (i = NPEVMIN; i < MAXPEV; i++)
    {
      PEv[i].filled = FALSE;
      PEv[i].EarlyGate = FALSE;
      PEv[i].val = 0;
      PEv[i].vsn = 0;
      PEv[i].ch = 0;
      bzero(PEv[i].name, 32);
    };
  };

  /* 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, "nevents")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &nEvents);
      CheckNoArgs(nret, 2, str);
      printf("will sort a max of %i events\n", nEvents);
      fflush(stdout);

    }
    /*----------------------------*/
    /*----------------------------*/
    /* here we include the common */
    /* options we have with GTSort */
    /*----------------------------*/

#include "GSSortReadChat_comoptions.h"

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

    else if ((p = strstr(str, "evtime0")) != NULL)
    {
      nret = sscanf(str, "%s %f", str1, &f1);
      EvTime0=(double)f1;
      printf("will start event time from %f seconds\n", (float) EvTime0);
      fflush(stdout);

    } else if ((p = strstr(str, "input")) != NULL)
      if (!InputSourceOnCommandLine)
      {
	nret = sscanf(str, "%s %s %s", str1, str2, str3);
	CheckNoArgs(nret, 3, str);
	if (strcmp("tape", str2) == 0)
	{
	  printf("will take input from tape\n");
	  set_gsII_indrive(str3);
	  InputSrc = TAPE;
	  fflush(stdout);
	} else if (strcmp("disk", str2) == 0)
	{
	  printf("will take input from disk\n");
	  set_gsII_indrive(str3);
	  InputSrc = DISK;
	  fflush(stdout);
	} else if (strcmp("net", str2) == 0)
	{
	  printf("will take input from UDP sender\n");
	  GSudpPort = atoi(str3);
	  SetupUdpReceiver(GSudpPort);
	  printf("__using port %i\n", GSudpPort);
	  InputSrc = NET;
	  fflush(stdout);
	} else
	{
	  printf("unknown input option: %s\n", str2);
	  printf("aborting\n");
	  fflush(stdout);
	  exit(0);
	};
      } else
      {
	printf("input source specified on commandline, ignoring chatfile specification\n");
      }



    else if ((p = strstr(str, "bgorfoffset")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &BGORFOffset);
      CheckNoArgs(nret, 2, str);
      printf("will will add %i to tbgo-tac2\n", BGORFOffset);
      fflush(stdout);

    } else if ((p = strstr(str, "esidemin")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &esidemin);
      CheckNoArgs(nret, 2, str);
      printf("will require %i minimum eside signal\n", esidemin);
      fflush(stdout);

    } else if ((p = strstr(str, "ebgomin")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &ebgomin);
      CheckNoArgs(nret, 2, str);
      printf("will require %i minimum BGO signal\n", ebgomin);
      fflush(stdout);

    } else if ((p = strstr(str, "setrecordver")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &RecordVer);
      CheckNoArgs(nret, 2, str);
      printf("will overwrite GS RecordVer with %i\n", RecordVer);
      SetRecordVer_tape(RecordVer);
      SetRecordVer_net(RecordVer);
      SetRecordVer_disk(RecordVer);
      fflush(stdout);

    } else if ((p = strstr(str, "tac2subge")) != NULL)
    {
      nret = sscanf(str, "%s %f", str1, &tac2FactorGE);
      CheckNoArgs(nret, 2, str);
      subtac2GE = TRUE;
      printf("will subtract TAC2 from ge times\n");
      printf("__using factor %f on tac2 first\n", tac2FactorGE);
      fflush(stdout);

    } else if ((p = strstr(str, "geehicalibration")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      RdGeCalFile(str2, ehiGeOffset, ehiGeGain);
      ehiDoGainCor = 1;
      fflush(stdout);

    } else if ((p = strstr(str, "geesidecalibration")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      RdGeCalFile(str2, esideGeOffset, esideGeGain);
      esideDoGainCor = 1;
      fflush(stdout);

    } else if ((p = strstr(str, "geelocalibration")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      RdGeCalFile(str2, eloGeOffset, eloGeGain);
      eloDoGainCor = 1;
      fflush(stdout);

    } else if ((p = strstr(str, "getalignment")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      RdOffFile(str2, tGEoffset);
      fflush(stdout);

    } else if ((p = strstr(str, "bgotalignment")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);

      RdOffFile(str2, tBGOoffset);
      fflush(stdout);

    } else if ((p = strstr(str, "tac2subbgo")) != NULL)
    {
      nret = sscanf(str, "%s %f", str1, &tac2FactorBGO);
      CheckNoArgs(nret, 2, str);
      subtac2BGO = TRUE;
      printf("will subtract TAC2 from BGO times\n");
      printf("__using factor %f on tac2 first\n", tac2FactorBGO);
      fflush(stdout);

    } else if ((p = strstr(str, "undotac2sub")) != NULL)
    {
      nret = sscanf(str, "%s %f %f", str1, &undotac2sub_factor, &undotac2sub_offset);
      CheckNoArgs(nret, 3, str);
      undotac2sub = TRUE;
      printf("will undo TAC2 subtraction on ge times\n");
      printf("__by adding back tac2*%f - %f\n", undotac2sub_factor, undotac2sub_offset);
      fflush(stdout);

    } else if ((p = strstr(str, "echo")) != NULL)
    {
      echo = TRUE;
      printf("will echo chatscript instructions\n");

    } else if ((p = strstr(str, "mkgetemats")) != NULL)
    {
      nret = sscanf(str, "%s %i %f %f %i %f %f", str1,
		    &GeETmats_e_nn, &f1, &f2,
		    &GeETmats_t_nn, &f3, &f4);
      CheckNoArgs(nret, 7, str);
      GeETmats_e_lo = (double) f1;
      GeETmats_e_hi = (double) f2;
      GeETmats_t_lo = (double) f3;
      GeETmats_t_hi = (double) f4;
      GeETmats = TRUE;
      printf("will generate germanium T-vs-E matrices\n");
      printf("__E: %i channels, from %f to %f\n", GeETmats_e_nn, GeETmats_e_lo, GeETmats_e_hi);
      printf("__T: %i channels, from %f to %f\n", GeETmats_t_nn, GeETmats_t_lo, GeETmats_t_hi);
      d1 = NGE * GeETmats_e_nn * GeETmats_t_nn * sizeof(float);
      d1 /= 1000000;
      printf("will use %9.2f MBytes + overhead for matrices\n", d1);

    } else if ((p = strstr(str, "lecroy")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      SetFERAVSN(str2, FERATYPE0);
      i1 = FERATYPE0;
      printf("setting type %i (LeCroy) for VSNs: %s\n", i1, str2);
      fflush(stdout);

    } else if ((p = strstr(str, "silena")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      SetFERAVSN(str2, FERATYPE1);
      i1 = FERATYPE1;
      printf("setting type %i (Silena) for VSNs: %s\n", i1, str2);
      fflush(stdout);

    } else if ((p = strstr(str, "export:phillips")) != NULL)
    {
      nret = sscanf(str, "%s %s %i", str1, str2, &i2);
      CheckNoArgs(nret, 3, str);
      SetFERAVSN(str2, FERATYPE3);
      SetExportModNWords(str2, i2);
      i1 = FERATYPE3;
      printf("setting type %i (Export module) for VSNs: %s\n", i1, str2);
      printf("__expect %i FERAs to be serviced by this export module\n", i2);
      fflush(stdout);

    } else if ((p = strstr(str, "export:4208TDC")) != NULL)
    {
      nret = sscanf(str, "%s %s %i", str1, str2, &i2);
      CheckNoArgs(nret, 3, str);
      SetFERAVSN(str2, FERATYPE6);
      SetlongRangeTDCNWords(str2, i2);
      i1 = FERATYPE6;
      printf("setting type %i (long range TDC) for VSNs: %s\n", i1, str2);
      printf("__expect %i TDCs to be serviced by this export module\n", i2);
      fflush(stdout);

    } else if ((p = strstr(str, "ortec")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      SetFERAVSN(str2, FERATYPE2);
      i1 = FERATYPE2;
      printf("setting type %i (Ortec) for VSNs: %s\n", i1, str2);
      fflush(stdout);

    } else if ((p = strstr(str, "nferadebug")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &i1);
      CheckNoArgs(nret, 2, str);
      SetnFeraDebug(i1);
      printf("Will print first %i fera interpretation debug statements\n", i1);
      fflush(stdout);

    } else if ((p = strstr(str, "translateVSN")) != NULL)
    {
      nret = sscanf(str, "%s %i %i", str1, &i1, &i2);
      CheckNoArgs(nret, 3, str);
      SetFERAvsntable(i1, i2);
      printf("Will translate raw VSN %i to %i:\n", i1, i2);
      fflush(stdout);

    } else if ((p = strstr(str, "feraclock1")) != NULL)
    {
      nret = sscanf(str, "%s %s %f %f", str1, str2, &f1, &f2);
      CheckNoArgs(nret, 4, str);
      SetFERAVSN(str2, FERATYPE4);
      i1 = FERATYPE4;
      printf("setting type %i (clock1) for VSNs: %s\n", i1, str2);
      SetClockPar(1, f1, f2);
      fflush(stdout);

    } else if ((p = strstr(str, "feraclock2")) != NULL)
    {
      nret = sscanf(str, "%s %s %f %f", str1, str2, &f1, &f2);
      CheckNoArgs(nret, 4, str);
      SetFERAVSN(str2, FERATYPE5);
      i1 = FERATYPE5;
      printf("setting type %i (clock2) for VSNs: %s\n", i1, str2);
      SetClockPar(2, f1, f2);
      fflush(stdout);

    } else if ((p = strstr(str, "countpevlevcheck")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &CountPevLevCheck);
      CheckNoArgs(nret, 2, str);
      printf("will revaluate counting pevs %i times\n", CountPevLevCheck);
      fflush(stdout);

    } else if ((p = strstr(str, "printpevraw")) != NULL)
    {
      PrintRawPEv = TRUE;
      printf("will print pev's before conditions are applied\n", Beta);
      fflush(stdout);

    } else if ((p = strstr(str, "overridempevcheck")) != NULL)
    {
      mpevOverride = TRUE;
      printf("will override mpev operation sanity checks\n");
      fflush(stdout);

    } else if ((p = strstr(str, "ferapev")) != NULL)
    {
      st = nret = sscanf(str, "%s %i %i %s %i %f %f", str1, &i2, &i3, str2, &i4, &f1, &f2);
      CheckNoArgs(nret, 4, str);
      PEv[LenPEv].vsn = i2;
      PEv[LenPEv].ch = i3;
      bcopy(str2, PEv[LenPEv].name, 32);
      printf("_ assign fera ");
      printf("VSN=%i, CH=%i ", PEv[LenPEv].vsn, PEv[LenPEv].ch);
      printf("to name: <%s>[#%i]; ", PEv[LenPEv].name, LenPEv);
      PEv[LenPEv].counting = FALSE;
      PEv[LenPEv].ratespno = FALSE;
      switch (st)
      {
      case 4:
	PEv[LenPEv].offset = 0;
	PEv[LenPEv].gain = 1;
	PEv[LenPEv].thresh = 0;
	PEv[LenPEv].calib = FALSE;
	printf("no calibration\n");
	break;
      case 5:
	PEv[LenPEv].offset = f1;
	PEv[LenPEv].gain = f2;
	PEv[LenPEv].thresh = i4;
	PEv[LenPEv].calib = TRUE;
	printf("calib: %f +%f*ch\n", (float) PEv[LenPEv].offset, (float) PEv[LenPEv].gain);
	break;
      case 7:
	PEv[LenPEv].offset = f1;
	PEv[LenPEv].gain = f2;
	PEv[LenPEv].thresh = i4;
	PEv[LenPEv].calib = TRUE;
	printf("calib: %f +%f*ch\n", (float) PEv[LenPEv].offset, (float) PEv[LenPEv].gain);
	break;
      default:
	printf("arguments in <%s> not valid\n", str);
	printf("please fix and try again\n");
	exit(-1);
      };

      /* check that this fera channel hasn't been assigned already */

      for (i = 0; i < LenPEv; i++)
	if (PEv[i].vsn == PEv[LenPEv].vsn)
	  if (PEv[i].ch == PEv[LenPEv].ch)
	  {
	    printf("this FERA channel has already been assigned\n");
	    printf("please fix problem and try again\n\n");
	    exit(-1);
	  };

      LenPEv++;

    } else if ((p = strstr(str, "mkpev")) != NULL)
    {
      st = nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      bcopy(str2, PEv[LenPEv].name, 32);
      PEv[LenPEv].counting = FALSE;
      PEv[LenPEv].ratespno = FALSE;
      printf("associated pev: %i with name: %s\n", LenPEv, PEv[LenPEv].name);
      LenPEv++;
    } else if ((p = strstr(str, "setspecial")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, str2);
      CheckNoArgs(nret, 2, str);
      SetSpecial(str2);
      printf("did SetSpecial(\"%s\")\n", str2);
    } else if ((p = strstr(str, "mktimemask")) != NULL)
    {
      nret = sscanf(str, "%s %s %s", str1, TMask[NTimeMasks].name, TMask[NTimeMasks].strrange);
      CheckNoArgs(nret, 3, str);
      str_decomp(TMask[NTimeMasks].strrange, L14BITS, TMask[NTimeMasks].lot);
      printf("created time mask #%i: <%s> [%s]\n",
	     NTimeMasks,
	     TMask[NTimeMasks].name,
	     TMask[NTimeMasks].strrange);
      fflush(stdout);
      NTimeMasks++;
      if (NTimeMasks >= MAX2DS)
      {
	i1 = MAX2DS;
	printf("too many timemasks, max= %i\n", i1);
	printf("abort\n");
	exit(1);
      };

    } else if ((p = strstr(str, "firstfile")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &FirstFile);
      CheckNoArgs(nret, 2, str);
      printf("will position tape at file %i\n", FirstFile);
      fflush(stdout);

    } else if ((p = strstr(str, "nposwarn")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &i1);
      CheckNoArgs(nret, 2, str);
      SetNPosWarn(i1);

    } else if ((p = strstr(str, "binbgo")) != NULL)
    {
      BinBGO = TRUE;
      printf("will bin BGO spectra\n");

    } else if ((p = strstr(str, "nfiles")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &Nfiles);
      CheckNoArgs(nret, 2, str);
      printf("will sort %i tape files\n", Nfiles);
      fflush(stdout);

    } else if ((p = strstr(str, "writerawfera")) != NULL)
    {
      WriteRawFera = 1;
      printf("will dump raw fera data\n");
      fflush(stdout);

    } else if ((p = strstr(str, "ignoreexternal")) != NULL)
    {
      DoNotInterpExt = TRUE;
      printf("will not interpret external fera data\n");
      fflush(stdout);

    } else if ((p = strstr(str, "anglefile")) != NULL)
    {
      nret = sscanf(str, "%s %s", str1, AngleFileName);
      CheckNoArgs(nret, 2, str);
      printf("will read GS angles from file: %s\n", AngleFileName);
      HaveAngleFileName = 1;
      fflush(stdout);

    } else if ((p = strstr(str, "ggmatrix")) != NULL)
    {
      nret = sscanf(str, "%s %s %i %f %f", str1, str2,
		    &i1, &f1, &f2);
      CheckNoArgs(nret, 5, str);
      bcopy(str2, GG[nGamxGam].name, 32);

      GG[nGamxGam].len = i1;
      GG[nGamxGam].lo = (double) f1;
      GG[nGamxGam].hi = (double) f2;

      printf("will make gamxgam matrix with range: %f to %f; ",
	     GG[nGamxGam].lo, GG[nGamxGam].hi);
      printf("with %i channels\n", GG[nGamxGam].len);
      printf("__name: %s\n", GG[nGamxGam].name);
      GG[nGamxGam].prompt_cond = -1;	/* mark unset */
      GG[nGamxGam].dtbggm_cond = -1;	/* mark unset */

      nGamxGam++;
      if (nGamxGam >= MAX2DS)
      {
	i1 = MAX2DS;
	printf("you have reached the max number\n");
	printf("or 2D matrices: %i\n", i1);
	printf("about...\n");
	printf("\n");
	exit(1);
      };
      fflush(stdout);

    } else if ((p = strstr(str, "devmatrix")) != NULL)
    {
      nret = sscanf(str, "%s %s %s %s %i %i %i", str1,
		    str2, str3, str4, &i1, &i2, &i3);
      CheckNoArgs(nret, 7, str);
      PType = FindPEvMatNo(str2, &DevMatrix[nDevMatrix].dlo);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      PType = FindPEvMatNo(str3, &DevMatrix[nDevMatrix].dhi);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      bcopy(str4, DevMatrix[nDevMatrix].name, 32);
      DevMatrix[nDevMatrix].nn = i1;
      DevMatrix[nDevMatrix].ylo = i2;
      DevMatrix[nDevMatrix].yhi = i3;
      printf("devmatrix x from %s(pev#=%i) ", str2, DevMatrix[nDevMatrix].dlo);
      printf("to %s(pev#=%i)\n", str3, DevMatrix[nDevMatrix].dhi);
      printf("__name: <%s>\n", DevMatrix[nDevMatrix].name);
      printf("__y: %i channels ", DevMatrix[nDevMatrix].nn);
      printf("from %i to %i\n", DevMatrix[nDevMatrix].ylo, DevMatrix[nDevMatrix].yhi);
      nDevMatrix++;
      fflush(stdout);

    } else if ((p = strstr(str, "binpev1d")) != NULL)
    {
      nret = sscanf(str, "%s %s %i %f %f", str1,
		    str2, &i1, &f1, &f2);
      CheckNoArgs(nret, 5, str);

      PB1D[NPEv1DBIN].len = i1;
      PB1D[NPEv1DBIN].lo = (double) f1;
      PB1D[NPEv1DBIN].hi = (double) f2;

      /* find the pev number */

      PType = FindPEvMatNo(str2, &PB1D[NPEv1DBIN].pev);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };

      printf("will bin pev: %i ", PB1D[NPEv1DBIN].pev);
      printf("(%s) ", PEv[PB1D[NPEv1DBIN].pev].name);
      printf("in spectrum [%s]\n", PEv[PB1D[NPEv1DBIN].pev].name);
      printf("__in %i channels from ", PB1D[NPEv1DBIN].len);
      printf(" %f to %f\n", PB1D[NPEv1DBIN].lo, PB1D[NPEv1DBIN].hi);
      NPEv1DBIN++;
      fflush(stdout);

    } else if ((p = strstr(str, "ratespectrum")) != NULL)
    {
      nret = sscanf(str, "%s %s %i %f %f", str1,
		    str2, &i1, &f1, &f2);
      CheckNoArgs(nret, 5, str);

      bcopy(str2, ratesp[nratesp].name, 32);
      ratesp[nratesp].len = i1;
      ratesp[nratesp].lo = (double) f1;
      ratesp[nratesp].hi = (double) f2;
      ratesp[nratesp].incrmt = (double) i1 / (f2 - f1);

      /* find the pev number */

      PType = FindPEvMatNo(str2, &ratesp[nratesp].pev);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      PEv[ratesp[nratesp].pev].ratespno = nratesp;

      printf("will make rate spectrum %s_rate(%i)\n", ratesp[nratesp].name, nratesp);
      nratesp++;
      fflush(stdout);

    } else if ((p = strstr(str, "firstevent")) != NULL)
    {
	  nret = sscanf(str, "%s %i", str1, &i1);
	  CheckNoArgs(nret, 2, str);
	  firstEvent = i1;
	  printf("will skip %i events\n", firstEvent);
	  fflush(stdout);

    } else if ((p = strstr(str, "binpev2d")) != NULL)
    {
      nret = sscanf(str, "%s   %s %i %f %f     %s %i %f %f", str1,
		    str2, &i1, &f1, &f2, str3, &i2, &f3, &f4);
      CheckNoArgs(nret, 9, str);
      PB2D[NPEv2DBIN].lenx = i1;
      PB2D[NPEv2DBIN].lox = (double) f1;
      PB2D[NPEv2DBIN].hix = (double) f2;
      PB2D[NPEv2DBIN].leny = i2;
      PB2D[NPEv2DBIN].loy = (double) f3;
      PB2D[NPEv2DBIN].hiy = (double) f4;
      PType = FindPEvMatNo(str2, &PB2D[NPEv2DBIN].pevx);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      PType = FindPEvMatNo(str3, &PB2D[NPEv2DBIN].pevy);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      NPEv2DBIN++;
      fflush(stdout);

    } else if ((p = strstr(str, "pevgam2d")) != NULL)
    {
      nret = sscanf(str, "%s   %s %s   %i %f %f   %i %f %f", str1,
	       PEvG[NGamPEv2DBin].name, str2, &i1, &f1, &f2, &i2, &f3, &f4);
      CheckNoArgs(nret, 9, str);

      PEvG[NGamPEv2DBin].pevlen = i1;
      PEvG[NGamPEv2DBin].pevlo = (double) f1;
      PEvG[NGamPEv2DBin].pevhi = (double) f2;
      PEvG[NGamPEv2DBin].gamlen = i2;
      PEvG[NGamPEv2DBin].gamlo = f3;
      PEvG[NGamPEv2DBin].gamhi = f4;
      PEvG[NGamPEv2DBin].prompt_cond = -1;
      PType = FindPEvMatNo(str2, &PEvG[NGamPEv2DBin].pev);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      printf("will bin gammas vs pev: %i <%s>; ", PEvG[NGamPEv2DBin].pev, PEv[PEvG[NGamPEv2DBin].pev].name);
      printf("__x: %i, %f-%f; ", PEvG[NGamPEv2DBin].pevlen, PEvG[NGamPEv2DBin].pevlo, PEvG[NGamPEv2DBin].pevhi);
      printf("__y: %i, %f-%f\n", PEvG[NGamPEv2DBin].gamlen, PEvG[NGamPEv2DBin].gamlo, PEvG[NGamPEv2DBin].gamhi);
      NGamPEv2DBin++;
      fflush(stdout);

    } else if ((p = strstr(str, "modpev")) != NULL)
    {
      nret = sscanf(str, "%s %s ", str1, str2);
      CheckNoArgs(nret, 2, str);

      if (strcmp("addconst", str2) == 0)
      {
	nret = sscanf(str, "%s %s %s %f ", str1, str2,
		      str3, &f2);
	CheckNoArgs(nret, 4, str);
	mpev[nmpev].type = 1;
	PType = FindPEvMatNo(str3, &mpev[nmpev].pev1);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	mpev[nmpev].factor = f2;
	mpev[nmpev].pev2 = FALSE;
	printf("will add constant %f to pev: %i\n", mpev[nmpev].factor, mpev[nmpev].pev1);
      } else if (strcmp("mulconst", str2) == 0)
      {
	nret = sscanf(str, "%s %s %s %f", str1, str2,
		      str3, &f2);
	CheckNoArgs(nret, 4, str);
	mpev[nmpev].type = 2;
	PType = FindPEvMatNo(str3, &mpev[nmpev].pev1);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	mpev[nmpev].factor = f2;
	mpev[nmpev].pev2 = FALSE;
	printf("will multiply pev: %i with \n", mpev[nmpev].pev1, mpev[nmpev].factor);
      } else if (strcmp("addpev", str2) == 0)
      {
	nret = sscanf(str, "%s %s %s %f %f %s", str1, str2,
		      str3, &f1, &f2, str4);
	CheckNoArgs(nret, 6, str);
	mpev[nmpev].type = 3;
	PType = FindPEvMatNo(str3, &mpev[nmpev].pev1);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	mpev[nmpev].offset = f1;
	mpev[nmpev].factor = f2;

	PType = FindPEvMatNo(str4, &mpev[nmpev].pev2);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	printf("will add pev: %i (with factor %f, offset %f) to pev: %i\n",
	       mpev[nmpev].pev2, mpev[nmpev].factor, mpev[nmpev].offset, mpev[nmpev].pev1);

	/* auto detect a counting pev */

	PEv[mpev[nmpev].pev1].counting |= PEv[mpev[nmpev].pev2].counting;

      } else if (strcmp("mulpev", str2) == 0)
      {
	nret = sscanf(str, "%s %s %s %f %s", str1, str2,
		      str3, &f2, str4);
	CheckNoArgs(nret, 5, str);
	mpev[nmpev].type = 4;
	PType = FindPEvMatNo(str3, &mpev[nmpev].pev1);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	mpev[nmpev].factor = f2;
	PType = FindPEvMatNo(str4, &mpev[nmpev].pev2);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	printf("will multiply pev: %i (with factor %f) on pev: %i\n",
	       mpev[nmpev].pev2, mpev[nmpev].factor, mpev[nmpev].pev1);
      } else if (strcmp("divpev", str2) == 0)
      {
	nret = sscanf(str, "%s %s %s %f %s", str1, str2,
		      str3, &f2, str4);
	CheckNoArgs(nret, 5, str);
	mpev[nmpev].type = 5;
	PType = FindPEvMatNo(str3, &mpev[nmpev].pev1);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	mpev[nmpev].factor = f2;
	PType = FindPEvMatNo(str4, &mpev[nmpev].pev2);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	printf("will divide pev: %i (with factor %f) with pev: %i\n",
	       mpev[nmpev].pev1, mpev[nmpev].factor, mpev[nmpev].pev2);
      } else if (strcmp("countcond", str2) == 0)
      {
	nret = sscanf(str, "%s %s %s in %s", str1, str2, str3, str4);	/* redo */
	CheckNoArgs(nret, 4, str);
	mpev[nmpev].type = 6;
	FindCondNo(str3, &mpev[nmpev].cond);
	PType = FindPEvMatNo(str4, &mpev[nmpev].pev1);
	if (PType != PEVTYPE)
	{
	  printf("wrong PEv type\n");
	  exit(-1);
	};
	printf("will count condition %i (%s) ",
	       mpev[nmpev].cond, PEvCond[mpev[nmpev].cond].name);
	printf("in pev: %i (%s)\n",
	       mpev[nmpev].pev1, PEv[mpev[nmpev].pev1].name);
	PEv[mpev[nmpev].pev1].counting = TRUE;
      } else
      {
	printf("second argument <%s> not understood\n", str2);
	printf("in chatscript line %i argument: %s\n", nn, str);
	exit(1);
      };

      nmpev++;
      if (nmpev >= MAXPEV)
      {
	printf("too many mpev statements...\n");
	exit(1);
      };

    } else if ((p = strstr(str, "cleanmultlim")) != NULL)
    {
      nret = sscanf(str, "%s %i %i", str1, &CleanMultLo, &CleanMultHi);
      CheckNoArgs(nret, 3, str);
      printf("will require ge clean multiplicity to be between %i and %i\n", CleanMultLo, CleanMultHi);
      fflush(stdout);

    } else if ((p = strstr(str, "dirtymultlim")) != NULL)
    {
      nret = sscanf(str, "%s %i %i", str1, &DirtyMultLo, &DirtyMultHi);
      CheckNoArgs(nret, 3, str);
      printf("will require ge dirty multiplicity to be between %i and %i\n", DirtyMultLo, DirtyMultHi);
      fflush(stdout);

    } else if ((p = strstr(str, "totalmultlim")) != NULL)
    {
      nret = sscanf(str, "%s %i %i", str1, &TotalMultLo, &TotalMultHi);
      printf("will require ge total multiplicity to be between %i and %i\n", TotalMultLo, TotalMultHi);
      fflush(stdout);

    } else if ((p = strstr(str, "isomerdata")) != NULL)
    {
      nret = sscanf(str, "%s mod %i segments %i %i %i %i %i %i %i", str1, &i1, &j1, &j2, &j3, &j4, &j5, &j6, &j7);
      CheckNoArgs(nret, 9, str);
      printf("will enable module %i for isomer tag data\n", i1);

      /* note: in principle we could aso have */
      /* and ignore/dontcare value here!      */

      /* set info */

      IsomerData = TRUE;

      setIsomerIDs(i1);;
      IsomerIDValid[i1][1] = j1;
      IsomerIDValid[i1][2] = j2;
      IsomerIDValid[i1][3] = j3;
      IsomerIDValid[i1][4] = j4;
      IsomerIDValid[i1][5] = j5;
      IsomerIDValid[i1][6] = j6;
      IsomerIDValid[i1][7] = j7;

      printf("_: ");
      for (i = 1; i <= 7; i++)
	switch (IsomerIDValid[i1][i])
	{
	case 0:
	  printf("[s%i,off]", i);
	  break;
	case 1:
	  printf("[s%i, on]", i);
	  break;
	case -1:
	  printf("[s%i,ign]", i);
	  break;
	default:
	  printf("IsomerIDValid[%i][%i]=%i not valid\n", i1, i, IsomerIDValid[i1][i]);
	  printf("you must specify 0,1,-1 for off and on and ignore\n");
	  exit(1);
	  break;
	}
      printf("\n");
      fflush(stdout);

    } else if ((p = strstr(str, "recalcgssume")) != NULL)
    {
      nret = sscanf(str, "%s %i %i %i %i", str1, &HKgelo, &HKgehi, &HKBGOlo, &HKBGOhi);
      CheckNoArgs(nret, 5, str);
      mkHK = TRUE;
      printf("will calculate HK in GSSort/overwrite EFF calculation\n");
      printf("__see pevs: \"sumge\" and \"sumbgo\"\n");
      printf("__ge  time limits from %4i to %4i\n", HKgelo, HKgehi);
      printf("__BGO time limits from %4i to %4i\n", HKBGOlo, HKBGOhi);
      fflush(stdout);

    } else if ((p = strstr(str, "nferaerrtoprint")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &NFERAErrToPrint);
      CheckNoArgs(nret, 2, str);
      printf("will print %i FERA errors at a time\n", NFERAErrToPrint);
      fflush(stdout);

    } else if ((p = strstr(str, "ignorehs")) != NULL)
    {
      ignoreHS = TRUE;
      printf("will ignore honeycomb suppression\n");

      fflush(stdout);

    } else if ((p = strstr(str, "gerfoffset")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &geRFOffset);
      CheckNoArgs(nret, 2, str);
      printf("will will add %i to tge-tac2\n", geRFOffset);
      fflush(stdout);

    } else if ((p = strstr(str, "hit2dmat")) != NULL)
    {
      nret = sscanf(str, "%s %s %s %s %s %s %i", str1,
		    str2, str3, str4, str5, str6, &i1);
      CheckNoArgs(nret, 7, str);
      PType = FindPEvMatNo(str2, &Hit2DMat[nHit2DMat].dlo1);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      PType = FindPEvMatNo(str3, &Hit2DMat[nHit2DMat].dhi1);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      PType = FindPEvMatNo(str4, &Hit2DMat[nHit2DMat].dlo2);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      PType = FindPEvMatNo(str5, &Hit2DMat[nHit2DMat].dhi2);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      bcopy(str6, Hit2DMat[nHit2DMat].name, 32);
      Hit2DMat[nHit2DMat].thresh = i1;
      fflush(stdout);

      /* make associated max energy pev */

      sprintf(str, "%s_emax", Hit2DMat[nHit2DMat].name);
      bcopy(str, PEv[LenPEv].name, 32);
      PEv[LenPEv].counting = FALSE;
      PEv[LenPEv].ratespno = FALSE;
      printf("__associated max E pev: #%i has name name: <%s>\n", LenPEv, PEv[LenPEv].name);
      Hit2DMat[nHit2DMat].emaxpev = LenPEv;	/* pev # to keep max energy
						 * in */
      LenPEv++;

      nHit2DMat++;

    } else if ((p = strstr(str, "compress")) != NULL)
    {
      nret = sscanf(str, "%s %i", str1, &ComPressLevel);
      CheckNoArgs(nret, 2, str);
      printf("will compress output file at level %i \n", ComPressLevel);
      fflush(stdout);
    } else if ((p = strstr(str, "pevcond1d")) != NULL)
    {
      nret = sscanf(str, "%s  %s  %s   %s", str1,
		    str2, PEvCond[NPEvCond].strrange, str3);
      CheckNoArgs(nret, 4, str);

      PEvCond[NPEvCond].type = 1;
      PEvCond[NPEvCond].filled = TRUE;
      str_decomp(PEvCond[NPEvCond].strrange, L14BITS, PEvCond[NPEvCond].lot);
      PType = FindPEvMatNo(str2, &PEvCond[NPEvCond].pev1);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type, you must apply to PEv\n");
	exit(-1);
      };
      PEvCond[NPEvCond].pev2 = 0;
      bzero(PEvCond[NPEvCond].fn, 72);
      bcopy(str3, PEvCond[NPEvCond].name, 32);
      PEvCond[NPEvCond].win = NULL;

      printf("made 1D condition # %i, <%s> , on pev %s, [%s]\n",
	 NPEvCond, PEvCond[NPEvCond].name, PEv[PEvCond[NPEvCond].pev1].name,
	     PEvCond[NPEvCond].strrange);

      NPEvCond++;
      if (NPEvCond >= MAXPEV)
      {
	printf("too many conditions\n");
	exit(1);
      };
      fflush(stdout);

    } else if ((p = strstr(str, "pevcond2d")) != NULL)
    {
      nret = sscanf(str, "%s  %s %s %s %s", str1,
		    str2, str3, str4, str5);
      CheckNoArgs(nret, 5, str);

      PType = FindPEvMatNo(str2, &PEvCond[NPEvCond].pev1);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      PType = FindPEvMatNo(str3, &PEvCond[NPEvCond].pev2);
      if (PType != PEVTYPE)
      {
	printf("wrong PEv type\n");
	exit(-1);
      };
      bcopy(str4, PEvCond[NPEvCond].fn, 72);
      bcopy(str5, PEvCond[NPEvCond].name, 32);

      PEvCond[NPEvCond].type = 2;
      PEvCond[NPEvCond].filled = TRUE;

      /* read in 2Dwindow */

      PEvCond[NPEvCond].win = NULL;
      PEvCond[NPEvCond].win = rd2dwin(PEvCond[NPEvCond].fn);

      printf("made 2D condition # %i, <%s> ", NPEvCond, PEvCond[NPEvCond].name);
      printf("on pev %s and %s\n", PEv[PEvCond[NPEvCond].pev1].name, PEv[PEvCond[NPEvCond].pev2].name);

      NPEvCond++;
      fflush(stdout);

    } else if ((p = strstr(str, "applycond")) != NULL)
    {
      nret = sscanf(str, "%s  %s  to  %s", str1,
		    str2, str3);
      CheckNoArgs(nret, 3, str);

      /* identify pseudo event vector or matrix to apply condition to */

      PType = FindPEvMatNo(str3, &ApplyCond[nApplyCond].pev);

      /* mark application target */

      ApplyCond[nApplyCond].pevtype = PType;

      /* identify this conditions number */

      i1 = 0;
      for (k = 0; k < NPEvCond; k++)
	if (strcmp(PEvCond[k].name, str2) == 0)
	{
	  i1++;
	  ApplyCond[nApplyCond].cond = k;
	};
      if (i1 == 1)
	printf("identified condition <%s> with condition number %i\n",
	       str2, ApplyCond[nApplyCond].cond);
      else if (i1 == 0)
      {
	printf("condition <%s> not found\n", str2);
	exit(0);
      } else if (i1 > 1)
      {
	printf("condition <%s> not unique\n", str2);
	exit(0);
      };

      /* confirm */

      printf("will apply condition #%i <%s> to ", ApplyCond[nApplyCond].cond, str2);
      if (PType == PEVTYPE)
	printf("pev %i\n", ApplyCond[nApplyCond].pev);
      else if (PType == GGMATTYPE)
	printf("gg matrix #%i named <%s>\n",
	     ApplyCond[nApplyCond].pev, GG[ApplyCond[nApplyCond].pev].name);
      else if (PType == GAMPEV)
	printf("ggampev matrix #%i\n", ApplyCond[nApplyCond].pev);

      ApplyCond[nApplyCond].inherited = FALSE;
      nApplyCond++;
      fflush(stdout);

    } else if ((p = strstr(str, "applyrftimemask")) != NULL)
    {
      nret = sscanf(str, "%s  %s  to  %s", str1,
		    str2, str3);
      CheckNoArgs(nret, 3, str);

      /* identify time mask */

      FindTimeMaskNo(str2, &i2);

      /* identify pseudo event vector or matrix to apply condition to */

      PType = FindPEvMatNo(str3, &i1);

      /* this better be a GG matrix! */

      if (!(PType == GGMATTYPE || PType == GAMPEV))
      {
	printf("you can only apply to a gamgam or gampev matrix\n");
	exit(1);
      };

      if (PType == GGMATTYPE)
      {

	/* mark RF time mask */

	if (GG[i1].prompt_cond != -1)
	{
	  printf("a prompt condition has already been apllied..\n");
	  exit(-1);
	};
	GG[i1].prompt_cond = i2;

	printf("apply (prompt) RF time mask #%i: <%s> [%s] to matrix #%i: <%s>\n",
	       GG[i1].prompt_cond, TMask[i2].name, TMask[i2].strrange, i1, GG[i1].name);
      };

      if (PType == GAMPEV)
      {

	/* mark RF time mask */

	if (PEvG[i1].prompt_cond != -1)
	{
	  printf("a prompt condition has already been apllied..\n");
	  exit(-1);
	};
	PEvG[i1].prompt_cond = i2;

	printf("apply (prompt) RF time mask #%i: <%s> [%s] to matrix #%i: <%s>\n",
	       PEvG[i1].prompt_cond, TMask[i2].name, TMask[i2].strrange, i1, GG[i1].name);
      };

    } else if ((p = strstr(str, "mkdssdmap")) != NULL)
    {
      nret = sscanf(str, "%s  %s  %s", str1,
		    str2, str3);
      CheckNoArgs(nret, 3, str);

      printf("will make DSSD FERA PeVs using mapfile: <%s>\n", str3);
      printf("__having names as %s###\n", str2);

      /* open DSSD map file */

      if ((fp1 = fopen(str3, "r")) == NULL)
      {
	printf("could not find <%s>\n", str3);
	printf("please fix problem and try again\n");
	exit(1);
      };
      printf("<%s> open\n", str3);

      /* read through the map file */

      pc1 = fgets(str4, STRLEN, fp1);
      while (pc1 != NULL)
      {

	/* extract map and calibration */

	nret = sscanf(str4, "%i %i %i %i %f %f %s %i %i", &i1, &i2, &i3, &i4, &f1, &f2, str5, &i5, &i6);

	if (nret == 6 || nret == 9)
	{

	  /* assign pev */

	  PEv[LenPEv].vsn = i1;
	  PEv[LenPEv].ch = i2;
	  sprintf(PEv[LenPEv].name, "%s%3.3i", str2, i3);
	  PEv[LenPEv].counting = FALSE;
	  PEv[LenPEv].ratespno = FALSE;
	  PEv[LenPEv].thresh = i4;
	  PEv[LenPEv].offset = (double) f1;
	  PEv[LenPEv].gain = (double) f2;
	  PEv[LenPEv].calib = TRUE;
	  PEv[LenPEv].EarlyGate = FALSE;

	  for (i = 0; i < LenPEv; i++)
	    if (PEv[i].vsn == PEv[LenPEv].vsn)
	      if (PEv[i].ch == PEv[LenPEv].ch)
	      {
		printf("this FERA channel has already been assigned\n");
		printf("please fix problem and try again\n\n");
		exit(-1);
	      };
	  printf("_DSSD FERA ");
	  printf("VSN=%i/CH=%i: ", PEv[LenPEv].vsn, PEv[LenPEv].ch);
	  printf("<%s>[#%i]; ", PEv[LenPEv].name, LenPEv);
	  printf("t=%4i,o=%6.2f g=%9.6f ", PEv[LenPEv].thresh,
		 (double) PEv[LenPEv].offset, (double) PEv[LenPEv].gain);

	  if (nret == 6)
	    printf(" -- no early gating\n");
	  else
	    printf("\n");

	};

	if (nret == 9)
	{

	  PEv[LenPEv].EarlyGate = TRUE;
	  PEv[LenPEv].EarlyGate_lo = (double) i5;
	  PEv[LenPEv].EarlyGate_hi = (double) i6;

	  /* identify the gating pev number */
	  /* and check it exists and is unique */

	  i1 = 0;
	  for (k = 0; k < LenPEv; k++)
	    if (strcmp(PEv[k].name, str5) == 0)
	    {
	      i1++;
	      PEv[LenPEv].EarlyGate_pev = k;
	    };
	  if (i1 == 0)
	  {
	    printf("gating pev %s does not exist\n");
	    exit(1);
	  } else if (i1 > 1)
	  {
	    printf("gating pev %s is not unique\n");
	    exit(1);
	  } else
	  {
	    printf("__ primary gated with pev= %s[#%i] ", str5, PEv[LenPEv].EarlyGate_pev);
	    printf("lo= %f to ", (float) PEv[LenPEv].EarlyGate_lo);
	    printf("hi= %f\n", (float) PEv[LenPEv].EarlyGate_hi);
	  };

	};

	if (!(nret == 6 || nret == 9))
	{
	  printf("problem with line\n");
	  printf("__ %s \n", str4);
	  printf("__ in file %s\n", str3);
	  exit(1);
	};

	/* increment the pev vector length */

	LenPEv++;

	/* get next line from map file */

	pc1 = fgets(str4, STRLEN, fp1);

      };

      /* close DSSD map file */

      fclose(fp1);


    } else if ((p = strstr(str, "applydttimemask")) != NULL)
    {
      nret = sscanf(str, "%s  %s  to  %s", str1,
		    str2, str3);
      CheckNoArgs(nret, 3, str);

      /* identify time mask */

      FindTimeMaskNo(str2, &i2);

      /* identify pseudo event vector or matrix to apply condition to */

      PType = FindPEvMatNo(str3, &i1);

      /* this better be a GG matrix! */

      if (PType != GGMATTYPE)
      {
	printf("you can only apply to a gamgam matrix\n");
	exit(1);
      };

      /* mark dt time mask */

      if (GG[i1].dtbggm_cond != -1)
      {
	printf("a delta time for gam condition has already been apllied..\n");
	exit(-1);
      };
      GG[i1].dtbggm_cond = i2;

      printf("apply gamma dt time mask #%i: <%s> [%s] to matrix #%i: <%s>\n",
	     GG[i1].dtbggm_cond, TMask[i2].name, TMask[i2].strrange, i1, GG[i1].name);

    }
    /* allow user to add chat script entries here */

#include "UserChat.h"

    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);
    /* rmEndComment(str, STRLEN); */

  };

  /* 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);

}

#if(0)
/*---------------------------------------------------------------------------*/

void
rmEndComment(char *str, int len)

{

  /* declarations */

  int             i;

  i = 0;
#if(0)
  while (str[i] != 0)
    printf("%3i:%i,%c\n", i, str[i], str[i++]);
#endif

};
#endif


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

void
WriteStatistics(unsigned int CurEvNo, int nMultOK, int *GSSortError, double runTime)
{

  /* declarations */

  int             i, i1, dt;
  int             EvProc;	/* events processed */
  float           r1;
  void            PrintFERAErrors(int);
  void            ReportBufPosErrors(void);
  int             err_print(int);
  int             hh,mm,ss;
  double          d1;

  /* tell all */

  printf("runtime: %f seconds, or ", runTime);
  d1=runTime;
  hh=(int)(d1/3600);
  d1-=(hh*3600);
  mm=(int)(d1/60);
  d1-=(mm*60);
  ss=d1;
  printf("%2ih%2im%2is\n", hh,mm,ss);
  
  printf("event counter %10i (total calls)\n", CurEvNo);
  EvProc=CurEvNo - firstEvent;
  printf("found         %10i good events\n", GSSortError[SUCCESS]); 
  printf("event rate: %9.3f ev/sec ", (float)GSSortError[SUCCESS]/(float)runTime);
  printf("[based on usec clock]\n");
  printf("processed     %10i buffers\n", GSSortError[GSIIERRNEWBUF]);
  r1=(float)GSSortError[SUCCESS]/(float)GSSortError[GSIIERRNEWBUF];
  printf("average number of events per buffer: %6.3f\n", r1);
  r1=(float)(DATASIZE-11)/r1;
  printf("roughly %9.3f words per event\n", r1);
  printf("processed %i KBytes\n", GSSortError[GSIIERRNEWBUF] * 16);
  fflush(stdout);
  dt = time(NULL) - tStart;
  printf("in %i seconds, integrated rates:\n", dt);
  printf("__%9.2f events/sec; ", (float) EvProc / (float) dt);
  r1 = 1000 * (float) dt / (float) GSSortError[GSIIERRNEWBUF];
  printf("%9.2f msec/ev; ", r1);
  r1 = (float) GSSortError[GSIIERRNEWBUF] / (float) dt;
  r1 *= 16.384;
  printf("%9.2f Kbytes/sec\n", r1);

  /* error/info statistics */

  for (i = 0; i < NERR; i++)
    if (GSSortError[i] > 0)
    {
      r1 = (float) GSSortError[i] * 100.0 / (float) EvProc;
      err_print(i);
      printf("%6.2f%% (%10i)\n", r1, GSSortError[i]);
#if(0)
      switch (i)
      {
      case GOODEV:
	printf("GOODEV\n");
	break;
      case GSIIERRNEWBUF:
	printf("GSIIERRNEWBUF\n");
	break;
      case GSIIERRFILEHEAD:
	printf("GSIIERRFILEHEAD\n");
	break;
      case GSIIGETNEWBUFF:
	printf("GSIIGETNEWBUFF\n");
	break;
      case GSIITAPEHEAD:
	printf("GSIITAPEHEAD\n");
	break;
      case GSIIMULTOOLOW:
	printf("GSIIMULTOOLOW\n");
	break;
      case GSIIERRSTART:
	printf("GSIIERRSTART\n");
	break;
      case GSIIERRSPACE:
	printf("GSIIERRSPACE\n");
	break;
      case GSIIERRMAXSANE:
	printf("GSIIERRMAXSANE\n");
	break;
      case GSIIERRBUFLEFT:
	printf("GSIIERRBUFLEFT\n");
	break;
      case GSIIERREBSIZE:
	printf("GSIIERREBSIZE\n");
	break;
      case FERAERRNOHEAD:
	printf("FERAERRNOHEAD\n");
	break;
      case FERAERRNEWHEAD:
	printf("FERAERRNEWHEAD\n");
	break;
      case GSIIERRCONT:
	printf("GSIIERRCONT\n");
	break;
      case GET_1_OF_N_RECVFROM_ERROR:
	printf("GET_1_OF_N_RECVFROM_ERROR\n");
	break;
      case GET_REST_OF_N_RECVFROM_ERROR:
	printf("GET_REST_OF_N_RECVFROM_ERROR\n");
	break;
      case EVENT_BUFFER_OVERFLOW_ERROR:
	printf("EVENT_BUFFER_OVERFLOW_ERROR\n");
	break;
      case UNEXPECTED_FRAGMENT_ERROR:
	printf("UNEXPECTED_FRAGMENT_ERROR (fragment drop)\n");
	break;
      case EMPTY_FRAGMENT_ERROR:
	printf("EMPTY_FRAGMENT_ERROR\n");
	break;
      case CHECKSUM_FAILURE_ERROR:
	printf("CHECKSUM_FAILURE_ERROR\n");
	break;
      default:
	printf("unknown error: %i\n", i);
	break;
      }
#endif
    };
  r1 = 100.0 * (float) nMultOK / (float) EvProc;
  printf("%6.2f%% passed the multiplicity requirements\n", r1);

  /* honeycomp suppression statistics */

  if (NoCleanWithoutHS > 0)
  {
    i1 = NoCleanWithoutHS + NoCleanWithHS;
    r1 = 100.0 * (float) NoCleanWithHS / (float) i1;
    printf("honeycomb suppression: %7.3f%% ", r1);
    r1 = 100.0 * (float) NoCleanWithHSmarked / (float) i1;
    printf("[specially marked: %7.3f%%]\n", r1);
  };

  /* simple statistics */

  i1 = NoCleanGE + NoDirtyGE + NoCleanBGO;
  r1 = 100.0 * (float) NoCleanGE / (float) i1;
  printf("NoCleanGE: %7.3f%% ", r1);
  r1 = 100.0 * (float) NoDirtyGE / (float) i1;
  printf("NoDirtyGE: %7.3f%% ", r1);
  r1 = 100.0 * (float) NoCleanBGO / (float) i1;
  printf("NoCleanBGO: %7.3f%%\n", r1);

  PrintFERAErrors(EvProc);

  ReportBufPosErrors();

#include "UserStat.h"

  /* done */

};

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

TCutG
* rd2dwin(Char_t * winname)

{

  /* declarations */

  char            str[STRLEN];
  TCutG          *mycutg;

  TFile          *f = new TFile(winname, "read");
  mycutg = (TCutG *) f->Get(winname);

  if (mycutg != NULL)
  {
    printf("2dwin read from file:\n");
    sprintf(str, "ls -l %s", winname);
    gSystem->Exec(str);
    /* mycutg->Print(); */
    fflush(stdout);
  } else
  {
    mycutg = NULL;
    printf("could not read 2dwin file %s\n", winname);
    printf("TFile error number: %i\n", f->GetErrno());
    printf("abort\n\n");
    exit(-1);
  };

  /* done */

  f->Delete();
  f->Close();
  return (mycutg);

}


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

void
FindTimeMaskNo(char *str, int *tm)

{
  /* declarations */

  int             k, i2, type;

  /* hunt for unique name */

  i2 = 0;
  for (k = 0; k < NTimeMasks; k++)
  {
    if (strcmp(TMask[k].name, str) == 0)
    {
      i2++;
      *tm = k;
    };
  };

  /* make sure the time mask name is unique */

  if (i2 > 1)
  {
    printf("time mask name %s is not unique\n", str);
    exit(0);
  } else if (i2 == 0)
  {
    printf("could not identify the time mask <%s>\n", str);
    printf("\n");
    printf("about\n");
    exit(-1);
  };

  /* done */

  return;

}

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

void
FindCondNo(char *str, int *nn)

{
  /* declarations */

  int             k, i2;

  /* hunt for unique name */

  i2 = 0;
  for (k = 0; k < NPEvCond; k++)
  {
    if (strcmp(PEvCond[k].name, str) == 0)
    {
      i2++;
      *nn = k;
    };
  };

  /* check that we found what we were looking for */
  /* -- or crash! */

  if (i2 == 0)
  {
    printf("cannot find %s in the condition list\n", str);
    exit(0);
  } else if (i2 > 1)
  {
    printf("condition name %s is not unique\n", str);
    exit(0);
  };

  /* done */

}

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

void
ApplyAllConditions()

{

  /* declarations */

  int             i;

  /* apply all conditions (to pevs and gg matrices alike) */

  for (i = 0; i < nApplyCond; i++)
    if (ApplyCond[i].pevtype == PEVTYPE)
    {
      if (PEv[ApplyCond[i].pev].filled && PEvCond[ApplyCond[i].cond].ok == 0)
	PEv[ApplyCond[i].pev].filled = 0;
    } else if (ApplyCond[i].pevtype == GGMATTYPE)
    {
      if (GG[ApplyCond[i].pev].filled && PEvCond[ApplyCond[i].cond].ok == 0)
	GG[ApplyCond[i].pev].filled = 0;
    };

}

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

void
FindConditions()
{

  /* declarations */

  int             i, i1;

  /* evaluate all conditions */

  for (i = 0; i < NPEvCond; i++)
    if (PEvCond[i].filled)
      switch (PEvCond[i].type)
      {
      case 1:
	i1 = (int) (PEv[PEvCond[i].pev1].val + 0.5);
	/* if (PEv[PEvCond[i].pev1].val >= PEvCond[i].lo &&
	 * PEv[PEvCond[i].pev1].val <= PEvCond[i].hi) */
	if (PEvCond[i].lot[i1])
	{
	  if (!PEvCond[i].ok)
	    PEvCond[i].count++;
	  PEvCond[i].ok = 1;
	} else
	  PEvCond[i].ok = 0;
	break;
      case 2:
	if (PEvCond[i].win->IsInside(PEv[PEvCond[i].pev1].val, PEv[PEvCond[i].pev2].val))
	{
	  if (!PEvCond[i].ok)
	    PEvCond[i].count++;
	  PEvCond[i].ok = 1;
	} else
	  PEvCond[i].ok = 0;
	break;

      };

};

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

void
FindParentPEv(int i)
{

  /* recursively find the parents of   */
  /* all pevs and their dependencies   */
  /* and apply relevant conditions     */
  /* from the inheritance tree         */

  /* ...this is 1/2 spooky!            */

  int             j, k;

  for (j = 0; j < nmpev; j++)
    if (mpev[j].type == 3 && mpev[j].pev1 == i)
    {
      printf("__parent: {pev:%i,name:%s} ", mpev[j].pev2, PEv[mpev[j].pev2].name);

      /* find conditions on this parent */

      printf("cond: ");
      for (k = 0; k < nApplyCond; k++)
	if (ApplyCond[k].pev == mpev[j].pev2)
	{
	  printf("[%s]", PEvCond[ApplyCond[k].cond].name);

	  /* apply this condition to last sibling as well */

	  ApplyCond[nApplyCond].inherited = TRUE;
	  ApplyCond[nApplyCond].pev = i;
	  ApplyCond[nApplyCond].cond = ApplyCond[k].cond;
	  ApplyCond[nApplyCond].pevtype = ApplyCond[k].pevtype;
	  nApplyCond++;

	}
      printf("\n");

      /* trace all the way back */

      if (PEv[mpev[j].pev2].vsn == 0)
	FindParentPEv(mpev[j].pev2);
    };

}

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

void
mpevOP_reg()
{

  /* do the mpev operations on the regular pevs */

  /* declarations */

  char            ok;
  int             i;
  double          d1;

  /* loop through the requested operations */

  for (i = 0; i < nmpev; i++)
    switch (mpev[i].type)
    {
    case 1:

      /* add constant */

      ok = PEv[mpev[i].pev1].filled;
      ok = ok || mpevOverride;
      if (ok)
      {
	PEv[mpev[i].pev1].filled = TRUE;
	PEv[mpev[i].pev1].val += mpev[i].factor;
      };
      break;

    case 2:

      /* multiply by constant */

      ok = PEv[mpev[i].pev1].filled;
      ok = ok || mpevOverride;
      if (ok)
      {
	PEv[mpev[i].pev1].filled = TRUE;
	PEv[mpev[i].pev1].val *= mpev[i].factor;
      };
      break;

    case 3:

      /* add another pev */

      ok = PEv[mpev[i].pev1].filled || PEv[mpev[i].pev2].filled;
      ok = ok || mpevOverride;
      ok = ok || PEv[mpev[i].pev1].counting;
      if (ok)
      {
	PEv[mpev[i].pev1].filled = TRUE;
	d1 = mpev[i].factor * PEv[mpev[i].pev2].val + mpev[i].offset;
	PEv[mpev[i].pev1].val += d1;
      };
      break;

    case 4:

      /* multiply by another pev */

      ok = PEv[mpev[i].pev1].filled && PEv[mpev[i].pev2].filled;
      ok = ok || mpevOverride;
      if (ok)
      {
	PEv[mpev[i].pev1].filled = TRUE;
	PEv[mpev[i].pev1].val *= (mpev[i].factor * PEv[mpev[i].pev2].val);
      };
      break;

    case 5:

      /* divide by another pev */


      ok = PEv[mpev[i].pev1].filled && PEv[mpev[i].pev2].filled;
      ok = ok || mpevOverride;
      if (ok)
      {
	PEv[mpev[i].pev1].filled = TRUE;
	PEv[mpev[i].pev1].val /= (mpev[i].factor * PEv[mpev[i].pev2].val);
      };
      break;

    };

};
/*---------------------------------------------------------------------------*/

void
mpevOP_special()
{

  /* do the mpev operations on the special zero pevs */

  /* declarations */

  char            ok;
  int             i;
  char            check[MAXPEV];
  double          d1;

  /* find all counting pevs that are currently zero */

  for (i = 0; i < LenPEv; i++)
    check[i] = FALSE;
  for (i = 0; i < LenPEv; i++)
    if (PEv[i].counting)
      if (PEv[i].val < 0.01)
	check[i] = TRUE;

#if(0)
  for (i = 0; i < LenPEv; i++)
    if (check[i])
    {
      printf("(%i,", i);
      printf("%s) ", PEv[i].name);
    };
  printf("\n");
#endif

  /* loop through the requested operations */

  for (i = 0; i < nmpev; i++)
    switch (mpev[i].type)
    {

    case 3:

      /* add pev */

      ok = check[mpev[i].pev1];
      if (ok)
      {
	PEv[mpev[i].pev1].filled = TRUE;
	d1 = mpev[i].factor * PEv[mpev[i].pev2].val + mpev[i].offset;
	PEv[mpev[i].pev1].val += d1;
      };
      break;

    default:
      break;

    };

};

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

void
CountConditions()
{

  /* declarations */

  int             i;

  for (i = 0; i < nmpev; i++)
    if (mpev[i].type == 6)
    {
      PEv[mpev[i].pev1].filled = TRUE;
      if (PEvCond[mpev[i].cond].ok)
	PEv[mpev[i].pev1].val++;
    };

}

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

void
MarkCountingPevs()

{

  /* find all pevs that are used */
  /* for counting and mark them  */

  /* all primary counting pevs already    */
  /* marked via 'modpev countcond'        */
  /* statements. Here we catch the ones   */
  /* made via modpev addpev operations    */

  /* declarations */

  int             i, ncp1, ncp2;

  /* first count all counting pevs */

  ncp1 = 0;
  for (i = 0; i < LenPEv; i++)
    if (PEv[i].counting)
      ncp1++;

  /* propagate counting flag */

  for (i = 0; i < nmpev; i++)
    if (mpev[i].type == 3)
      if (PEv[mpev[i].pev2].counting)
	PEv[mpev[i].pev1].counting = TRUE;

  /* count again */

  ncp2 = 0;
  for (i = 0; i < LenPEv; i++)
    if (PEv[i].counting)
      ncp1++;

  /* shall we look again ? (recursively) */
  /* (keep trying until we can't find anymore */

  if (ncp2 > ncp1)
    MarkCountingPevs();

}


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

void
PrintCondToApply()
{
  int             i;

  printf("\n");
  printf("primary+inhereted conditions to be applied:\n");
  printf("\n");
  for (i = 0; i < nApplyCond; i++)
  {
    printf("%3i: apply condition # %3i, <%s>, ",
	   i, ApplyCond[i].cond, PEvCond[ApplyCond[i].cond].name);

    /* list the object, which can be different kinds */

    if (ApplyCond[i].pevtype == PEVTYPE)
      printf("to pev # %3i named <%s>", ApplyCond[i].pev, PEv[ApplyCond[i].pev].name);
    else if (ApplyCond[i].pevtype == GGMATTYPE)
      printf("to gg matrix # %3i named <%s>", ApplyCond[i].pev, GG[ApplyCond[i].pev].name);
    else if (ApplyCond[i].pevtype == GAMPEV)
      printf("to gampev matrix # %3i named <%s>", ApplyCond[i].pev, PEvG[ApplyCond[i].pev].name);
    else
    {
      printf("unexpected ApplyCond[i].pevtype.... failing\n");
      exit(-1);
    }

    /* more info */

    if (ApplyCond[i].inherited)
      printf(", INHERITED, ");
    else
      printf(", PRIMARY  , ");

    /* redundant object type */

    if (ApplyCond[i].pevtype == PEVTYPE)
      printf("pev   type ");
    else if (ApplyCond[i].pevtype == GGMATTYPE)
      printf("ggmat type ");
    else if (ApplyCond[i].pevtype == GAMPEV)
      printf("gampev type ");
    printf("\n");

  };
  printf("\n");
};

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

int
FindPEvMatNo(char *str, int *num)

/* find an object number */

{
  /* declarations */

  int             k, i2, type;

  /* hunt for unique name */

  i2 = 0;


  for (k = 0; k < LenPEv; k++)
  {
    if (strcmp(PEv[k].name, str) == 0)
    {
      i2++;
      *num = k;
      type = PEVTYPE;
    };
  };

  for (k = 0; k < nGamxGam; k++)
  {
    if (strcmp(GG[k].name, str) == 0)
    {
      i2++;
      *num = k;
      type = GGMATTYPE;
    };
  };

  for (k = 0; k < NGamPEv2DBin; k++)
  {
    if (strcmp(PEvG[k].name, str) == 0)
    {
      i2++;
      *num = k;
      type = GAMPEV;
    };
  };

  /* make sure the name is unique */

  if (i2 > 1)
  {
    printf("object name <%s> is not unique\n", str);
    printf("abort\n");
    exit(0);
  };

  /* quit if we do not find the name */

  if (i2 == 0)
  {
    printf("no such object name: <%s>\n", str);
    printf("please fix problem and try again\n");
    exit(-1);
  };

  /* done */

  return (type);

}



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

void
DoPrimaryGating()
{

  /* declarations */

  int             i;

  /* loop through and primary gate */

  for (i = 0; i < LenPEv; i++)
    if (PEv[i].EarlyGate)
    {
      if (PEv[PEv[i].EarlyGate_pev].val > PEv[i].EarlyGate_hi)
	PEv[i].filled = FALSE;
      if (PEv[PEv[i].EarlyGate_pev].val < PEv[i].EarlyGate_lo)
	PEv[i].filled = FALSE;
    };

  /* done */

};

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

void
SetBeta()
{

  /* delarations */

  int             i;
  double          d1;

  /*-------------------------------------*/
  /* find Doppler and aberration factors */
  /*-------------------------------------*/

  DopCorFac[0] = 0;
  ACFac[0] = 0;
  for (i = 1; i < NGE; i++)
  {
    printf("det %3.3i-> ", i);
    d1 = angle[i] / 57.29577951;
    DopCorFac[i] = (1 - Beta * cos(d1)) / sqrt(1 - Beta * Beta);
    DopCorFacHires[i] = DopCorFac[i];
    printf("dop cor fac: %6.4f; ", DopCorFac[i]);
    ACFac[i] = DopCorFac[i] * DopCorFac[i];
    printf("aberration cor fac: %6.4f\n", ACFac[i]);

  };
  fflush(stdout);

  /*-------------------------------------------*/
  /* add hi res data compression via DopCorFac */
  /*-------------------------------------------*/

  for (i = 1; i < NGE; i++)
    DopCorFacHires[i] *= HiResDataMultFactor;
  printf("hires data compression set via DopCorFacHires[]\n");
  fflush(stdout);

  /* done */

}

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

int
decode_isomer(unsigned int *ext,
	      unsigned int *next,
	      unsigned int *xtype,
	      unsigned int *xvsn,
	      unsigned int *xch,
	      unsigned int *xdata,
	      unsigned short int tac2)

/* input */

/* extdata[];	isomer data in ext FERA array */

/* output */

/* next;	 # isomer modules found */
/* xtype[];	 LeCroy/Silena type array         */
/* xvsn[];	 VSN array                        */
/* xch[];	 channel array                    */
/* xdata[];	 data array                       */

{
  /* declarations */

  int             nn, i, j, seg[8], k, nseg, ngood, nbad, nignore, nfera, i1;
  static int      bincode[15] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512,
  1024, 2048, 4096, 8192, 16384};
  unsigned short int hitpat;
  unsigned int    id, tt, ee;

  /* interpret data */

  nn = ext[0] / 3;
  nfera = 0;
  *next = 0;

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

    j = (3 * i) + 1;		/* ext pos */

    /* find ID and hitpattern */

    id = ext[j] & 0x00ff;
    modsum[id]++;
#if(0)
    if (id > 7)
    {
      for (k = 0; k <= ext[0]; k++)
	printf("ext[%3i]= %i = 0x%4.4x\n", k, ext[k], ext[k]);
    };
#endif
    hitpat = (ext[j] & 0xfe00) >> 9;
    nseg = nbad = ngood = nignore = 0;
    for (k = 0; k < 8; k++)
      if ((hitpat & bincode[k]) == bincode[k])
      {
	seg[nseg++] = k + 1;	/* count as 1...7 */
	if (IsomerIDValid[id][k + 1] == 1)
	  ngood++;
	else if (IsomerIDValid[id][k + 1] == 0)
	  nbad++;
	else if (IsomerIDValid[id][k + 1] == -1)
	  nignore++;
	segsum[id][k + 1]++;
      };

    /* if good event, xfer data to fera arrays */
    /* what is good is not totally unique perhaps... */

    if (ngood > 0 && nbad == 0)
    {

      /* extract time and energy */

      if (subtac2BGO)
      {

	i1 = (int) (ext[++j] & M12BITS);
	i1 -= (int) (tac2 * tac2FactorBGO);
	i1 += BGORFOffset;
	if (i1 >= 0 && i1 < M12BITS)
	  tt = (unsigned short int) i1;
	else
	  tt = 0;

      } else
	tt = ext[++j] & M12BITS;

#if(0)
      if (NprintEvNo <= NumToPrint)
	if (tt < 0 || tt > 3000)
	{
	  printf("ext[j-1]=%i\n", ext[j - 1]);
	  printf("(int) (tac2 *tac2FactorBGO)=%i\n", (int) (tac2 * tac2FactorBGO));
	  printf("BGORFOffset=%i\n", BGORFOffset);
	  printf("i1=%i\n", i1);
	  printf("tt=%i\n", tt);
	};
#endif

      ee = ext[++j] & M12BITS;


      /* return as interpreted fera data */
      /* if time and energy is ok        */

      if (tt > 0 && ee > 0)
      {
	xtype[nfera] = FERATYPE5;
	xvsn[nfera] = id;
	xch[nfera] = 1;
	xdata[nfera++] = ee;
	xtype[nfera] = FERATYPE5;
	xvsn[nfera] = id;
	xch[nfera] = 2;
	xdata[nfera++] = tt;
      };

    };

    /* print interpretation */

    if (NprintEvNo <= NumToPrint)
    {
      printf("iso mod: %3.3i, nseg: %2.2i; ", id, nseg);
      printf("ngood=%i, nbad=%i, nignore=%i -> ", ngood, nbad, nignore);
      if (ngood > 0 && nbad == 0)
      {
	printf("t: %4.4i, ", tt);
	printf("(%4.4i) ", ext[j - 1]);
	printf("e: %4.4i; ", ee);
	printf("nfera=%i; ", nfera);
      };
      printf("\n");
      fflush(stdout);
    };



  };

  /* count succesful isomer tags */

  nallIsomerEv++;
  if (nfera > 0)
    ngoodIsomerEv++;

  /* done */

  NprintEvNo++;
  *next = nfera;
  return (SUCCESS);

}


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

void
wrIsomerStatistics(int ntot)

{

  /* This function is called at the end */
  /* of the sort to present the */
  /* statistics collected by the */
  /* decode_isomer function. */

  /* declarations */

  double          sum, d1;
  int             i, j;

  printf("enter \"wrIsomerStatistics\"\n");

  /* total statistics */

  printf("\n");
  d1 = 100 * (double) ngoodIsomerEv / (double) ntot;
  printf(" raw Isomer tag event fraction: %9.3f%%\n", d1);
  d1 = 100 * (double) ngoodIsomerEv / (double) ntot;
  printf("good Isomer tag event fraction: %9.3f%%\n", d1);

  /* module statistics */

  sum = 0;
  for (i = 0; i < NGE; i++)
    sum += modsum[i];

  printf("\n");
  printf("module statistics\n");
  printf("\n");
  for (i = 0; i < NGE; i++)
    if (modsum[i] > 0)
    {
      printf("isomer mod %3i: ", i);
      printf("counts: %10i, ", modsum[i]);
      printf(" %9.2f%%\n", 100 * (double) modsum[i] / sum);
    };

  printf("\n");
  printf("segment statistics\n");
  printf("\n");
  for (i = 0; i < NGE; i++)
    if (modsum[i] > 0)
    {
      sum = 0;
      for (j = 1; j <= 7; j++)
	sum += segsum[i][j];
      for (j = 1; j <= 7; j++)
      {
	printf("isomer mod %3i: ", i);
	printf("seg: %1i; ", j);
	printf("counts: %10i, ", segsum[i][j]);
	printf(" %9.2f%%; ", 100 * (double) segsum[i][j] / sum);
	if (IsomerIDValid[i][j] == 1)
	  printf(" <-isomer segment\n");
	else if (IsomerIDValid[i][j] == 0)
	  printf(" <-veto segment\n");
	else if (IsomerIDValid[i][j] == -1)
	  printf(" <-ignore segment\n");
      };
      printf("\n");
    };

  /* Tac2 subtraction info/warning */

  printf("\n");
  if (subtac2BGO)
    printf("Tac2 was subtracted from BGO times, factor: %f\n", tac2FactorBGO);
  else
  {
    printf("WARNING: tac2 was not subtracted from BGO times!\n");
    printf("_use: \"tac2subbgo factor\" in chat script\n");
  };

  printf("exit \"wrIsomerStatistics\"\n");

};

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

#include "GSGTComFunctions.h"