/*
* Copyright (C) 2002-2009, CompHEP Collaboration
* Author: Alexander Sherstnev 
* ------------------------------------------------------
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

#include "service2/include/syst.h"

#include "tag_reader.h"
#include "tag_writer.h"
#include "tag_parser.h"
#include "fill.h"
#include "validity.h"
#include "structures.h"
#include "compare.h"
#include "mix_format2.h"

int mix_format2 (const int nf, const char target[], const char names[], const int lenth, long the_seed)
{
  int jf, t, tt, nx;
  int nSubtot;
  int nProc;
  int *shift;
  long i,l;
  long nRectot;
  long maxnRec;
  long num;
  long **pos;
  long fpos;
  long *pos0;
  long *nRec;
  long *real_nRec;
  long *nLeft;
  double *ri, xrn;
  double sigmatot = 0.0;
  char buff[1024];
  char buff1[1024];
  char **map;
  char **filename;

  FILE** infile;
  FILE* outFile;

  process_ theProcessUP;
  tags **tbase;
  proc_pos *proc_position;

  srand48 (the_seed);

/****************************************************************/
/*  allocate the memory space for work dimensions */
  tbase    = malloc (nf * sizeof (tags *));
  pos0     = malloc (nf * sizeof (long));
  shift    = malloc (nf * sizeof (int));
  infile   = malloc (nf * sizeof(FILE*));
  filename = malloc (nf * sizeof(char*));

  for (i = 0; i < nf; ++i) {
    filename[i] = malloc ((lenth + 1) * sizeof(char));
    strncpy (filename[i], names + i * lenth, lenth);
    filename[i][lenth] = 0;
    trim (filename[i]);
  }

/****************************************************************/
/* check arguments - input event files                          */
  fprintf (stdout, "mix (info): files to mix and randomize:\n");
  for (i = 0; i < nf; ++i) {
    infile[i] = fopen (filename[i], "r");
    if (!infile[i]) {
      fprintf (stdout, "mix (error): file not found: %s\n", filename[i]);
      return 2;
    } else {
      fprintf (stdout, "%s\n", filename[i]);
    }
  }

  for (i = 0; i < nf; ++i) {
    tbase[i] = init_cap(100);
    cup_reader (infile[i], tbase[i]);
    syntax_validity (tbase[i]);
    structure_validity (tbase[i]);
    physics_validity (tbase[i]);
    pos0[i] = ftell (infile[i]) + 1;  /* position of the 1st event beginning */
  }

  fill_LH_structures (nf, tbase, &theProcessUP, names, lenth);
  compare_extra_info (nf, tbase, &theProcessUP);
  nSubtot = theProcessUP.proc_info.NprocRUP;
  proc_position = malloc (nSubtot * sizeof (proc_pos));
  compare_processes (nf, tbase, proc_position);

/****************************************************************/
/*   allocate the memory space for dimensions                   */
  nRec = malloc (nSubtot * sizeof (long));
  real_nRec = malloc (nSubtot * sizeof (long));
  nLeft = malloc (nSubtot * sizeof (long));
  ri = malloc (nSubtot * sizeof (double));
  pos = malloc (nSubtot * sizeof (long *));
  map = malloc (nSubtot * sizeof (char *));

  nRectot = 0;
  maxnRec = 0;
  for (i = 0; i < nSubtot; i++)
    real_nRec[i] = 0;

  t = 0;
  for (i = 0; i < nf; i++) {
    int nEvCount = 0;
    tt = get_tag (0, tbase[i], "total");
    if (tt != -1) {
      tt = get_ival (0, "Nproc", tbase[i]->tag[tt]);
    } else {
      fprintf (stderr, "mix (error): the file %s does not contain the total tag\n", filename[i]);
      exit (2);
    }

    fseek (infile[i], pos0[i], SEEK_SET);
    fprintf (stdout, "mix (info): reading file %s:\n", filename[i]);
    while (1) {
      if (!fgets (buff, 1000, infile[i]))
        break;
      sscanf (buff, "%i", &nProc);
      num = t + nProc - 1;
      if (num < nSubtot) {
        ++real_nRec[num];
      } else {
        fprintf (stderr, "mix (error): process number of event exceed max process number in file %s\n", filename[i]);
        exit (2);
      }
      if (0 == nEvCount % 20000) fprintf (stdout, "mix (info): %i events read\n", nEvCount);
      ++nEvCount;
    }
    fprintf (stdout, "mix (info): %i events read\n", nEvCount);
    t += tt;
  }

  sigmatot = 0.0;
  for (i = 0; i < nSubtot; i++) {
    nRec[i] = proc_position[i].Nevents;
    if (real_nRec[i] != nRec[i]) {
      fprintf (stdout, " mix (warning): real number of events in file %s (%li) does not coinside with header info (%li)\n",
               filename[i], real_nRec[i], nRec[i]);
      nRec[i] = real_nRec[i];
    }
    nLeft[i] = nRec[i];
    nRectot = nRectot + nRec[i];
    if (nRec[i] > maxnRec) {
      maxnRec = nRec[i];
    }
    pos[i] = malloc (nRec[i] * sizeof (long));
    map[i] = malloc (nRec[i] * sizeof (char));
    sigmatot += theProcessUP.proc_info.crossecUP[i];
  }

  t = 0;
  for (i = 0; i < nSubtot; i++) {
    real_nRec[i] = 0;
  }
  for (i = 0; i < nf; i++) {
    tt = get_tag (0, tbase[i], "total");
    if (tt != -1) {
      tt = get_ival (0, "Nproc", tbase[i]->tag[tt]);
    } else {
      fprintf (stderr, "mix (error): the file %s des not contain the total tag.\n", filename[i]);
      exit (2);
    }
    shift[i] = t;

    fseek (infile[i], pos0[i], SEEK_SET);
    while (1) {
      fpos = ftell (infile[i]);
      if (!fgets (buff, 1000, infile[i]))
        break;
      sscanf (buff, "%i", &nProc);
      num = t + nProc - 1;
      pos[num][real_nRec[num]] = fpos;
      map[num][real_nRec[num]] = 0;
      real_nRec[num]++;
    }
    t += tt;
  }

/******************************************************************/
/* evaluate the map of subprocesses weights in the interval [0,1] */
  xrn = 0;
  for (i = 0; i < nSubtot; i++) {
    ri[i] = xrn + theProcessUP.proc_info.crossecUP[i] / sigmatot;
    xrn = ri[i];
  }

/****************************************************************/
/*         mix events from the files                            */
  outFile = fopen (target, "w");

  write_cap (outFile, tbase, theProcessUP);

  for (i = 0; i < nRectot; i++)
  {
    xrn = drand48 ();
    jf = 0;
    while (xrn > ri[jf]) {
      jf++;
    }

    if (nLeft[jf] == 0) break;
    l = 0;
    while (map[jf][l]) {
      l++;
      if (l == nRec[jf]) {
        l = 0;
      }
    }
    fseek (infile[proc_position[jf].nfile], pos[jf][l], SEEK_SET);
    map[jf][l] = 1;
    fgets (buff, 1000, infile[proc_position[jf].nfile]);
    sscanf (buff, "%i:%[^\n]", &nx, buff1);
    fprintf (outFile, " %i:%s\n", (int) (shift[proc_position[jf].nfile] + nx), buff1);
    --(nLeft[jf]);
  }

  fclose (outFile);

  fprintf (stdout, "\nmix (info): final statistics:\n");
  for (i = 0; i < nf; ++i) {
    fprintf (stdout,  "mix (info): file %s: %li events used from %li events", filename[i], nRec[i] - nLeft[i], nRec[i]);
    if (0 == nLeft[i])  fprintf (stderr, ", file exhausted!\n");
    else fputs ("\n", stdout);
  }
  fprintf (stdout, "mix (info): %li mixed/randomized events written to %s\n", nRectot, target);

  final_write_cap (target, nLeft, i);
  free (pos);
  free (map);
  free (proc_position);
  free (tbase);
  free (pos0);
  free (shift);
  free (nRec);
  free (nLeft);
  free (ri);
  for (i = 0; i < nf; ++i) {
    free (filename[i]);
  }
  free (filename);

  return 1;
}