/*****************************************************
 *     generate_hb_test_events.c
 *
 * C Program to generate one file with Nev HB test events.
 * Assumes HB maps are produced by generate_hb_maps, so
 * it simply uses the same algorithm.
 * These two program must be kept syncrhonised, this is ugly,
 * but allows for simple fast code which shold be OK for simple
 * fast, though extensive, tests.
 * A mininal variation is allowed, as this program reads the
 * AM_MAP from disk, so allowing to change at will the road
 * assignement in generate_hb_maps.c as long as the SS list and
 * the hit->SS correspondance is the same.
 * This is done to avoid a lenghty inversion of the SS_MAP
 * This program needs to know the SuperStrip size, 
 * to decide if InternalOverlfow is generated and to figure out
 * the maximum number of SuperStrips per layer, so it reads it
 * form the same SS_size file as generate_hb_maps.c does.
 *
 * S.Belforte 16 Oct, 1999.
 *
 * revision history:
 * 18 Oct 99: VERSION 2: change arguments to:
 *                  RANDOM (0=false, 1=true)
 *                  Niter =1,2 only used for RANDOM=false to
 *                         divide systematic test in smaller files
 *                  Multiplicity = only used for RANMDON=false to
 *                                 set the #of hits for each SS, 
 *                                 which is always the same
 *  Version 2 frozen on 19 Oct 1999. No major mileston, but most
 *      of random stuff is there and hopefully without having
 *      messed up the simpler systematics checks.
 *      It is frozen to allow more free work to complete
 *      and debug random code for Version 3.
 *
 *****************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "/home/cdf_svt/lib/vmesvt/error_def.h"
#include "/home/cdf_svt/lib/vmesvt/tools.h"

/* GLOBAL DEFINITIONS AND VARIABLES */

const unsigned long EEflag = 0x400000; /* bit 22 */
const unsigned long EPflag = 0x200000; /* bit 21 */
const unsigned long IOVFflag = 0x2000; /* bit 13 */

/**************************************************************************/
/*** if following constans change, code has to change !!!!!!!!!!! *********/
const int Max_Iter = 2; /* number of iterations for full systematic check */
const int MaxExtraWords = 1; /* maximum extra words for hit, i.e. 2 total */
const int MaxNoisy =  2; /* maximum number of noisy hits per layer        */
const int MaxRoadsInEvent = 20; /* maximum number of roads in one event   */
/**************************************************************************/

int Random, iter, multiplicity;
unsigned long * ss_am_map;
unsigned long * hb_am_map;
int SSsize;

int get_SS_size();
void parse_arguments(int argc, char* argv[]);
int pick_number_of_events();
int pick_number_of_roads();
unsigned long find_SS  (unsigned long road, unsigned long layer);
int find_hits (unsigned long SS, unsigned long layer,
	       unsigned long * hits, int* n2out, unsigned long* IOVF);
int extra_word();
int pick_noise(int layer);
int kill_hit ();
int Parity(unsigned long w, unsigned long n);
int iran (int imin, int imax);
int rnpoissn (float mean);

/**************** CODE *********************/

int main (int argc, char *argv[])
{
 int N_SS;    /* number of SuperStrip per layer */
 int MaxHitsInLayer;
 int MaxRoadsInEvent;
 int evmin, evmax;
 int Nev, ev;
 int Nroad, roads_in_this_event;

 int hb_am_map_len;
 int i, nhit,nout,nroad;
 int hitparity, roadparity, outparity;
 int words, bytes;
 unsigned long road, layer, SSid, SS, hit, out;
 int hit_words;
 int n2out;
 unsigned long * hits_in_SS;
 unsigned long bunch_id;
 unsigned long L2buff;
 unsigned long L1bits;
 unsigned long IOVF, intoverflow;
 unsigned long data;
 unsigned long * Hits;
 unsigned long * Roads;
 unsigned long * Outs;
 FILE* ss_map;
 FILE* am_map;
 FILE* HITfile;
 FILE* ROADfile;
 FILE* OUTfile;

 N_SS = get_SS_size();
  
 parse_arguments(argc, argv);

 if (Random)
   {
     Nroad = 128 * 1024; /* 128k roads exists in AM Map */
     Nev = pick_number_of_events();
     evmin = 0;
     evmax = Nev;
   }
 else
   {
     Nev = N_SS/Max_Iter;
     evmin = Nev*(iter-1);
     evmax = Nev*(iter);
   }

 /* allocate memory, just use fixed arrays over and over */
 bytes = sizeof(unsigned long);
 MaxHitsInLayer = (SSsize + 1)*(1+MaxExtraWords) + MaxNoisy;
 words = MaxRoadsInEvent + 1;
 words = 10000;
 Roads = malloc (words*bytes);
 if (Roads == NULL) {printf("No memory for %d roads\n",words);exit(8);}
 words=MaxHitsInLayer;
 hits_in_SS = malloc(words*bytes);
 if (hits_in_SS == NULL) {printf("No memory for %d hits\n",words);exit(8);}
 words = MaxRoadsInEvent*8*MaxHitsInLayer + 1;
 words = 10000;
 Hits  = malloc (words*bytes);
 if (Hits == NULL) {printf("No memory for %d hits\n",words);exit(8);}
 words = MaxRoadsInEvent*(8*MaxHitsInLayer+1) + 1;
 words = 10000;
 Outs  = malloc (words*bytes);
 if (Outs == NULL) {printf("No memory for %d outs\n",words);exit(8);}

printf("Will make files with %d events \n",Nev);

/* read HB_AM_MAP */
 if (Random) {
   printf("reading HB_AM_MAP file .... ");
   fflush(NULL);
   hb_am_map_len=1024*1024; /* 1mega words */
   hb_am_map = (unsigned long *)malloc(hb_am_map_len * sizeof(unsigned long));
   hb_am_map_len=load_struct_hex("HB_AM_MAP",hb_am_map);
   printf(" DONE\n");
 }

/* open output files */
 if ( (HITfile = fopen("HIT_DATA", "w")) == NULL)
   { printf("cannot open file HIT_DATA for writing \n");  exit(0); }
 if ( (ROADfile = fopen("ROAD_DATA", "w")) == NULL)
   { printf("cannot open file ROAD_DATA for writing \n");  exit(0); }
 if ( (OUTfile = fopen("OUT_DATA", "w")) == NULL)
   { printf("cannot open file OUT_DATA for writing \n");  exit(0); }

/************* main loop *********************/

 for (ev=evmin; ev<evmax; ev++) {
   nhit  = 0;
   nout  = 0;
   nroad = 0;
   hitparity  = 0;
   roadparity = 0;
   outparity  = 0;
   intoverflow = 0;
   /* Start by generating ROAD_DATA file */
   /* pick a number of road to generate */
   if (!Random) 
     roads_in_this_event = 1;
   else
     roads_in_this_event = pick_number_of_roads();
   /* dimension data structure, maximum is SSsize +1 hits for
      each layer of each roads, plus maximum number of
      noisy hits plus 1 EE word */
   /* process each road: add to event record and find hits */
   for(i=0;i<roads_in_this_event;i++)
     {
       /* pick a road in the [0,Nroad) range, if not random test,
	  will just pick the ev-th road */
       if (!Random)
	 road = ev;
       else
	 road = iran(0,Nroad-1);
       if (!(ev%100)) printf ("processing road: %d\n", road);
       /* add road to road event record. Each road is a packet */
       Roads[nroad++] = EPflag | road;
       roadparity ^= Parity((EPflag|road),22);
       /* loop on all layers of this road */
       for(layer=0;layer<=7;layer++)
	 {
	   /* pick the SuperStrips from the AM map */
	   SS = find_SS (road, layer);
	   /* find the hits that produces this SuperStrip (could be >1 words)*/
	   hit_words = find_hits(SS, layer, hits_in_SS, &n2out, &IOVF);
	   intoverflow |= IOVF;
	   /* add hit to hit event record and the output record*/
	   for (i=0;i<hit_words;i++){
	     hit = hits_in_SS[i];
	     Hits[nhit++] = hit;
	     hitparity ^= Parity(hit,22);
	     if (i<n2out){
	       out = hit & 0x1FFFFF ; /*output hits are not packets*/
	       Outs[nout++] = out;
	       outparity ^= Parity(out, 22);
	     }
	   }
	 } /* end loop on layers */
       /* add end packet word with road number to out record */
       out = EPflag | road;
       Outs[nout++] = out;
       outparity ^= Parity(out, 22);
     } /* end loop on roads_in_this_event */
   /* compute End Event for Road and Hit and Out */
   /* L2buff and L1bits give 4 more bits for event identification ! */
   L2buff = ( (ev&0xC00) >>10 ) <<19;
   L1bits = ( (ev&0x300) >>8 )  <<17;
   bunch_id = ev & 0xFF;
   Hits[nhit]   = EEflag | EPflag | L2buff| L1bits | hitparity <<8 | bunch_id;
   Roads[nroad] = EEflag | EPflag | L2buff| L1bits | roadparity<<8 | bunch_id;
   Outs[nout]   = EEflag | EPflag | L2buff| L1bits | outparity <<8 | bunch_id;
   if (intoverflow) Outs[nout] |= IOVFflag ;
   /* write events to files and release memory */
   for (i=0;i<=nhit;i++)  fprintf(HITfile, "%06X\n", Hits[i]);
   for (i=0;i<=nroad;i++) fprintf(ROADfile,"%06X\n", Roads[i]);
   for (i=0;i<=nout;i++)  fprintf(OUTfile, "%06X\n", Outs[i]);
 } /* end of loop on events */
/* close files */
fclose (HITfile);
fclose (ROADfile);
fclose (OUTfile);

free (hits_in_SS);
free (Hits);
free (Roads);
free (Outs);

}  /* end of main */


/*************************** FUNCTIONS *******************************/


/*------------------------------------------------------------------*/
int get_SS_size()
{
  FILE* SSsize_file;
  int N_SS;

  if ( (SSsize_file = fopen("SS_size", "r")) == NULL)
    {
      printf("cannot open file SS_size \n");
      exit(0);
    }
  else
    {
      fscanf(SSsize_file,"%d", &SSsize);
      printf("file read, SSsize : %d \n", SSsize);
    }

  if (SSsize==2)
    N_SS = 4*1024;   /* 4k SS for each of 8 layers */
  else if (SSsize==4)
    N_SS = 2*1024;   /* 2k SS for each of 8 layers */
  else if (SSsize==8)
    N_SS =   1024;   /* 1k SS for each of 8 layers */
  else
    {printf(" Illegal SSsize specification %d. Abort\n", SSsize); exit(0);}

return N_SS;
}

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

void parse_arguments(int argc, char* argv[])
{
 if (argc<=1)
   {
     printf("Program Reads HB_AM_MAP and writes HIT_DATA, ROAD_DATA and OUT_DATA\n");
     printf("Syntax:\n\n");
     printf("%s <random_flag> <iteration_number> <multiplicity>\n\n",argv[0]);
     printf(" random_flag :  1 ==> fully random test, see code \n");
     printf("                0 ==> systematic check of all HLM locations \n");
     printf("\n");
     printf(" iteration_number :  ==> only used if random_flag=0 \n");
     printf("                     divide full systematic check in \n");
     printf("                     many parts to avoid overflowing files \n");
     printf("                     as of Oct 18 values are: 1-2 \n");
     printf("\n");
     printf(" multiplicity :  ==> only used if random_flag=0 \n");
     printf("                 a number in [0,SS_size+1] range sets \n");
     printf("                 the hit multiplicity in each SuperStrip \n");
     printf("\n");
     exit(0);
   }
  
 sscanf(argv[1], "%d", &Random);
 if (Random==1)
   printf("Will generate random events \n");
 else if (Random==0) {
   printf("Will generate one event for each Road \n");
   printf(" and systematically check all HLM locations \n"); }
 else
   { printf("Illegal Random_Flag value: %d. Abort \n",Random); exit(0);}

 if (!Random) {
   if (argc!=4) {
     printf("Missing required argument.\n");
     printf("Type \" %s \" for instructions\n",argv[0]);
     exit(0);
   }
   sscanf(argv[2], "%d", &iter);
   if ( (iter<=0) || (iter>Max_Iter) )
     { 
       printf("Illegal Iteration value: %d. Abort \n",iter);
       printf("Type \" %s \" for instructions\n",argv[0]);
       exit(0);
     }
   sscanf(argv[3], "%d", &multiplicity);
   if ( (multiplicity<0) || (multiplicity>(SSsize+1)) )
     { 
       printf("Illegal Multiplicity value: %d. Abort \n",multiplicity); 
       printf("Type \" %s \" for instructions\n",argv[0]);
       exit(0);
     }
 }
 return;
}

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

int pick_number_of_events()
{
  /* select a random numer of events for this file
     at present just a  Poisson with lower and upper limit */
const int MaxEv = 2000;
int nev;
#if 0
 do
   nev = rnpoissn((float)MaxEv/2.);
 while ( (nev==0) | (nev>MaxEv)) ;
#endif
nev = 1;
return nev;
}

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

int pick_number_of_roads()
{
  /* select a random numer of roads for this event.
     At present Poisson distribution with average 2.5
     This gives a peak at nroad=2 and about as many 1
     as 3. Also cut at MaxRoadsInEvent to be safe, less then 0.1%
     of the Poisson is above 10 anyhow  */

int nroad;
 do
   nroad = rnpoissn(2.5);
 while (nroad>MaxRoadsInEvent) ;
return nroad;
}

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

unsigned long find_SS (unsigned long road, unsigned long layer)
{
  /*
   */
unsigned long SSid;
unsigned long SS;
int am_map_add;

if (!Random)
  SS = road;
 else
   {
     am_map_add = road + (layer<<17);
     SSid = hb_am_map[am_map_add] & 0x1FFF; /* pick leftmost 13bits */
     if ((SSid&0x1)==0)  /* SSsize = 2 */
       SS = SSid >> 1;
     else if ((SSid&0x2)==0x1) /* SSsize = 4 */
       SS = SSid >> 2;
     else if ((SSid&0x2)==0x3) /* SSsize = 8 */
       SS = SSid >> 3;
     else { printf("Illegal SSid %x at AM_MAP add %x, wrong SSsize bits.",
		   SSid,road+layer<<13);
     printf("Abort\n"); exit(0); }
   }
return SS; 
}

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

int find_hits (unsigned long SS, unsigned long layer,
	       unsigned long * hits, int* n2out, unsigned long* IOVF)
{
  /* assume simple algorithm used in generate_hb_maps.c
   * So just left shit and pad with random numbers the
   * rightmost 4 bits that are dropped before addressing SS map
   * Note that layer has a different shift...
   * Also remember that each hit is a packet.
   */
  float hit_average;
  unsigned long right_byte;
  unsigned long SSdata;
  unsigned long SSid;
  unsigned long hit;
  int nwords, nhit, i;
  int noisy_hits;
  if(Random)
    {
      /* sometimes there are additional hits in the same SS
	 just use Poisson with average cut in (0,SSsize+1]
	 using average such to keep IOVF from occuring
	 in each event */
      switch(SSsize)
	{
	case 2: hit_average=0.5; break;
	case 4: hit_average=2.0; break;
	case 8: hit_average=5.0; break;
	}
      do
	nhit=rnpoissn(hit_average);
      while ((nhit==0)|(nhit>(SSsize+1))) ;
      /* find from SS Map one hit that belongs to this SS */
      /* for the time being do this as for non-random making
	 assumptions on how SS Map has been filled */
      SSdata = SS & 0x1FFF;
      /* create hits */
      nwords = 0;
      *IOVF = 0;
      for (i=0;i<nhit;i++)
	{
	  /* each hit has X% probability of missing */
	  if (kill_hit()) continue;
	  right_byte = iran(0,15);
	  hits[nwords++] = layer<<18 | SSdata<<4 | right_byte ;
	  /* sometime one hit is two words, fill with 16 random bits
	     and keep same layer, this just helps debugging... */
	  if (extra_word()) {
	    hits[nwords++] = layer<<18 | iran(0,64*1024) ;
	  }
	  /* mark last word of the packet */
	  hits[nwords-1] |= EPflag;
	} /* end loop on hits on the given SuperStrip */
      /* only as many words as teh HB can hold will go out */
      *n2out = nwords;
      if (*n2out>SSsize) *n2out=SSsize;
      /* tell caller if Internal Overflow will be generated */
      if (nwords>SSsize) *IOVF = 1;
      /* sometimes there are noisy hits in other SS of the same layer */
      noisy_hits = pick_noise(layer);
      for (i=0;i<noisy_hits;i++){
	do {
	  hit = layer<<18 | iran(0,64*1024); /* 16 random bits */
	  SSid = (hit & 0x1FFF) >> 4; /* SS this hit will point to */
	}
	while (SSid == SS); /* make sure noise is on other SS's */
	hit |= ( EPflag | (layer<<18) ) ; /* add layer and EP */
	hits[nwords++] = hit;
      }
    }
  else
    {         /* systematic check, just simple single word hits */
      nwords = multiplicity;
      SSdata = SS & 0x1FFF;
      for (i=0;i<nwords;i++)
	{
	  right_byte = iran(0,15);
	  hits[i] = EPflag | layer<<18 | SSdata<<4 | right_byte ;
	}
      if (nwords>SSsize) {
	*IOVF = 1;
	*n2out  = SSsize;
      }
      else {
	*IOVF = 0;
	*n2out  = nwords;
      }
    }
  return nwords; 
}

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

int kill_hit ()
{
  /* 10% of the hits are killed at present */
int test, kill;
kill=0;
test=iran(1,10);
if (test==1) kill=1;
return kill;
}

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

int extra_word ()
{
  /* 10% of hits have one extra word at present */
int test, extra;
extra =0;
test=iran(1,100);
if (test<=10) extra = 1;
return extra;
}

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

int pick_noise (int layer)
{
  /* decide how many extra hits on this layer, to be
     realistic this should be pretty large, to avoid
     being overflowed with irrelevant make a peak
     at 0 with a long flat tail. At present noise is the
     same on all layers */
int nnoisy, test;
test=iran(1,2);
if(test==1)
  { /* very little noise */
    do
      nnoisy=rnpoissn(0.2);
    while ( nnoisy>MaxNoisy ) ;
  }
else
  { /* long tail */
    nnoisy=iran(0,MaxNoisy);
  }
return nnoisy;
}


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

int Parity(unsigned long w, unsigned long n)
/*
 * Returns parity of the n LSb's of the word w. (Must be n>1). The
 * result is stored in the LSb of the return value. All other bits
 * are randomly set and should not be used from G. Punzi Parity.c
 */
{
  unsigned long   out = w;

  while (--n)  out ^= w >>= 1;
  out &= 0x1;
  return out; 
}

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

float ranf(void)
/* random in [0,1), 0 included, 1 excluded */
{
  float r;
  r = (float)rand() / ((float)(RAND_MAX)+1.);
  return r;
}

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

int iran (int imin, int imax) 
/*
 * Returns an integer random number
 * uniformly distributed in [imin,imax]
 * both imin and imax are a possible return value.
 */
{
  float r;      /* random in [0,1) */
  r = ranf();
  return imin+(int)(r*(float)(imax-imin+1));
}

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

int ulran (unsigned long ulmin, unsigned long ulmax) 
{
  float r; /* random in [0,1), 0 included, 1 excluded */
  r=1.;
  while (r==1.) r = ranf();
  return ulmin+(unsigned long)(r*(float)(ulmax-ulmin+1));
}


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

float fran() 
{
  float fmin=0.;
  float fmax=1.;
  float r;
  r=1.;
  while (r==1.) r=(float)(drand48()/1.);
  return (fmin+(r*(fmax-fmin)));
}

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

int rnpoissn(float mean)
/*
 * Returns interger distributed as Poisson with given mean
 * Translated from poiprb.F in CDF r_n area
 * by Bill AShmanskas on Oct 19 1999.
 */
{
  const float plim = 16; /* cut-off to switch to gaussian */
  const float TWOPI  = 2. * 3.141592653; /* just 2-pi */
  float a, p, r, rr, s, phi, x;
  int n;
  a = mean;
  n = 0;
  if (a>0) {
    if (a<=plim) {
      r = ranf();
      p = exp(-a);
      s = p;
      if (r>s) do {
        n++;
        p *= a/n;
        s += p;
      } while (s<r && p>1.0e-30);
    } else {
      rr = sqrt(-2*log(ranf()));
      phi = TWOPI*ranf();
      x = rr*cos(phi);
      n = a + x*sqrt(a) + 0.5;
      if (n<0) n = 0;
    }
  }
  return n;
}