/* 
* Copyright (C) 2001-2009, CompHEP Collaboration
* Copyright (C) 1997, Alexander Pukhov 
* ------------------------------------------------------
*/
#include "service2/include/chep_limits.h"
#include "service2/include/tptcmac.h"
#include "service2/include/read_func.h"
#include "service2/include/syst.h"
#include "chep_crt/include/edittab.h"
#include "chep_crt/include/crt_util.h"
#include "chep_crt/include/file_scr.h"
#include "plot/include/plot.h"
#include "out_ext.h"

#include "rd_num.h"
#include "subproc.h"
#include "phys_val.h"
#include "histogram.h"

table histTab = { "*** Table ***", "Distributions",
  "  Parameter  |> Min bound <|> Max bound <|> Rest Frame <|", NULL
};

typedef struct histRec
{
  struct histRec *next;
  linelist mother;
  long nPoints;
  char key;
  char plist[10];
  char restfr[10];
  double hMin, hMax;
  double f[300];
  double ff[300];
}
histRec;

static histRec *histPtr = NULL;

int
clearHists (void)
{
  histRec *hists = histPtr;
  int i;
  if (histPtr == NULL)
    return 0;
  while (hists)
    {
      for (i = 0; i < 300; i++)
        {
          hists->f[i] = 0;
          hists->ff[i] = 0;
        }
      hists->nPoints = 0;
      hists = hists->next;
    }
  return 1;
}



void
fillHists (double w)
{
  histRec *hists = histPtr;
  int i;
  double z;
  while (hists)
    {
      z = calcPhysVal (hists->key, hists->plist, hists->restfr);

      if (z > hists->hMin && z < hists->hMax)
        {
          i = 300 * (z - hists->hMin) / (hists->hMax - hists->hMin);
          hists->f[i] += w;
          hists->ff[i] += w * w;
        }
      hists->nPoints++;
      hists = hists->next;
    }
}


int
correctHistList (void)
{
  int lineNum = 0;
  int i = 0, n, k;
  char keyChar;
  double rMin;
  double rMax;
  char fieldName[50];

  midstr histStr;
  midstr minStr;
  midstr maxStr;
  midstr restStr;
  linelist ln = histTab.strings;
  histRec *hptr = histPtr;

  while (hptr != NULL)
    {
      hptr->mother = NULL;
      hptr = hptr->next;
    }

  while (ln != NULL)
    {
      char lv[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
      char restnumbers[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
      histStr[0] = 0;
      minStr[0] = 0;
      histStr[0] = 0;
      restStr[0] = 0;

      sscanf (ln->line, "%[^|]%*c%[^|]%*c%[^|]%*c%[^\n]%*c", histStr, minStr,
              maxStr, restStr);
      lineNum++;

      trim (minStr);
      trim (histStr);
      trim (maxStr);
      trim (restStr);

      strcpy (fieldName, "Wrong field 'Min. bound'");
      if (!minStr[0] || calcExpression (minStr, rd_num, &rMin))
        goto errorExit;

      strcpy (fieldName, "Wrong field 'Max bound'");
      if (!maxStr[0] || calcExpression (maxStr, rd_num, &rMax))
        goto errorExit;

      strcpy (fieldName, "Wrong field 'Parameter'");
      if (!checkPhysVal (histStr, &keyChar, lv))
        goto errorExit;

      if (rMax < rMin)
        {
          double tmp = rMax;
          rMax = rMin;
          rMin = tmp;
        }

/* New restframe */
      strcpy (fieldName, "Wrong field 'Rest frame'");
      if (restStr[0] == '|')
        {
          restnumbers[0] = 0;
        }
      else
        {
          i = 0;
          while (restStr[i] && restStr[i] != ' ')
            {
              n = restStr[i] - '0';

              if (n <= 0 || n > nin_ + nout_)
                goto errorExit;

              for (k = 0; k < i; k++)
                if (restnumbers[k] == n)
                  goto errorExit;

              restnumbers[i] = n;
              i++;
            }

          if (i < 10)
            restnumbers[i] = 0;

          while (i < 10)
            {
              if (restStr[i] && (restStr[i] != ' '))
                goto errorExit;
              i++;
            }
        }

      hptr = histPtr;
      while (hptr != NULL)
        {
          if ((hptr->key == keyChar) && (hptr->hMin == rMin)
              && (hptr->hMax == rMax) && (!strcmp (lv, hptr->plist))
              && (!strcmp (restnumbers, hptr->restfr)))
            {
              if (hptr->mother)
                {
                  strcpy (fieldName, " Duplicate record");
                  goto errorExit;
                }
              else
                {
                  hptr->mother = ln;
                  goto cont;
                }
            }
          else
            hptr = hptr->next;
        }
      hptr = malloc (sizeof (histRec));
      hptr->next = histPtr;
      hptr->mother = ln;
      histPtr = hptr;
      hptr->key = keyChar;
      strcpy (hptr->plist, lv);

      strcpy (hptr->restfr, restnumbers);

      hptr->hMin = rMin;
      hptr->hMax = rMax;
      for (i = 0; i < 300; i++)
        {
          hptr->f[i] = 0;
          hptr->ff[i] = 0;
        }
      hptr->nPoints = 0;
cont:
      ln = ln->next;
    }

  hptr = histPtr;
  histPtr = NULL;
  while (hptr != NULL)
    {
      histRec *hptr_ = hptr;
      hptr = hptr->next;
      if (!hptr_->mother)
        free (hptr_);
      else
        {
          hptr_->next = histPtr;
          histPtr = hptr_;
        }
    }

  return 0;

errorExit:
  sprintf (errorText, " Error in  line %d .\n" "%s.", lineNum, fieldName);
  warnanykey (2, 10, errorText);
  return 1;
}


static void
writeDistributions (FILE * iprt)
{
  histRec *hists = histPtr;
  int i;
  linelist rec;

  for (rec = histTab.strings; rec; rec = rec->next)
    {
      for (hists = histPtr; hists; hists = hists->next)
        if (hists->mother == rec)
          {
            fprintf (iprt, " nPoints=%ld key=%c", hists->nPoints, hists->key);
            for (i = 0; i < 10 && hists->plist[i]; i++)
              fprintf (iprt, "%d", hists->plist[i]);
            fprintf (iprt, " hMin=%-12E hMax=%-12E", hists->hMin,
                     hists->hMax);
            fprintf (iprt, "\n  f: ");
            for (i = 0; i < 300; i++)
              fprintf (iprt, " %-12E", hists->f[i]);
            fprintf (iprt, "\n ff: ");
            for (i = 0; i < 300; i++)
              fprintf (iprt, " %-12E", hists->ff[i]);
          }
      fprintf (iprt, "\n");
    }
}


static void
readDistributions (FILE * iprt)
{
  int i;
  linelist rec;
  for (rec = histTab.strings; rec; rec = rec->next)
    {
      histRec *hist = malloc (sizeof (histRec));

      hist->next = histPtr;
      hist->mother = rec;
      histPtr = hist;

      fscanf (iprt, " nPoints=%ld key=%c", &(hist->nPoints), &(hist->key));
      for (i = 0; i < 10; i++)
        {
          fscanf (iprt, "%c", hist->plist + i);
          if (hist->plist[i] == ' ')
            {
              hist->plist[i] = 0;
              break;
            }
          else
            hist->plist[i] -= '0';
        }
      fscanf (iprt, "hMin=%lf hMax=%lf\n", &(hist->hMin), &(hist->hMax));
      fscanf (iprt, "  f: ");
      for (i = 0; i < 300; i++)
        fscanf (iprt, " %lf", hist->f + i);
      fscanf (iprt, " ff: ");
      for (i = 0; i < 300; i++)
        fscanf (iprt, " %lf", hist->ff + i);
    }
}


int
WriteHistograms (FILE * nchan)
{
  fprintf (nchan, "\n");
  writetable0 (&histTab, nchan);
  writeDistributions (nchan);
  return 0;
}

int
ReadHistograms (FILE * nchan)
{
  fscanf (nchan, "\n");
  readtable0 (&histTab, nchan);
  readDistributions (nchan);
  return 0;
}

static int
nBinMenu (void)
{

  static int kBinMenu = 3;
  char strmen[] =
    "\015"
    " 300         "
    " 150         "
    " 100         "
    "  75         "
    "  60         "
    "  50         "
    "  30         "
    "  25         "
    "  20         "
    "  15         "
    "  12         "
    "  10         "
    "  6          "
    "  5          " "  4          " "  3          " "  2          ";

  void *pscr = NULL;

  int n;
  if (!kBinMenu)
    kBinMenu = 3;
  for (;;)
    {
      menu1 (54, 14, "number of bins", strmen, "", &pscr, &kBinMenu);
      if (kBinMenu)
        {
          sscanf (strmen + 1 + strmen[0] * (kBinMenu - 1), "%d", &n);
          put_text (&pscr);
          return n;
        }
      return 0;
    }
}


void
showHist (void)
{
  char histStr[STRSIZ], minStr[STRSIZ], maxStr[STRSIZ], restStr[STRSIZ];;
  linelist ln = histTab.strings;
  char *menutxt;
  void *pscr = NULL;
  int mode = 0;

  int npos = 0;
  int ii;
  int width, width2;
  char restnumbers[10] = { ' ' };

  while (ln)
    {
      npos++;
      ln = ln->next;
    }
  if (!npos)
    return;

  sscanf (histTab.format, "%[^|]%*c%[^|]%*c%[^|]%*c%[^\n]%*c", histStr,
          minStr, maxStr, restStr);

  width = strlen (histStr);
  width2 = strlen (restStr);

  menutxt = malloc (2 + (width + width2 + 3) * npos);
  menutxt[0] = width + width2 + 3;
  menutxt[1] = 0;

  ln = histTab.strings;
  while (ln)
    {
      sscanf (ln->line, "%[^|]%*c%[^|]%*c%[^|]%*c%[^\n]%*c", histStr, minStr,
              maxStr, restStr);
      strcat (menutxt, histStr);
      if ((restStr[0] == '|') || (restStr[0] == ' ') || (restStr[0] == '\n'))
        strcat (menutxt, "   ");
      else
        strcat (menutxt, "rf ");
      strcat (menutxt, restStr);
      ln = ln->next;
    }

  for (;;)
    {
      menu1 (54, 10, "", menutxt, "", &pscr, &mode);
      switch (mode)
        {
        case 0:
          free (menutxt);
          return;
        default:
          {
            histRec *hist = histPtr;
            int nBin;
            ln = histTab.strings;
            for (npos = 1; npos < mode; npos++)
              ln = ln->next;
            for (; hist && hist->mother != ln; hist = hist->next)
              {;
              }
            if (hist)
              {
                char xname[200], yname[80], units[80];
                if (hist->nPoints == 0)
                  messanykey (10, 10, "Distibution is empty");
                else
                  while ((nBin = nBinMenu ()))
                    {
                      double f[300], df[300], coeff;
                      int i;
                      coeff =
                        nBin / (hist->nPoints * (hist->hMax - hist->hMin));
                      for (i = 0; i < nBin; i++)
                        {
                          int k;
                          f[i] = 0;
                          df[i] = 0;
                          for (k = 0; k < 300 / nBin; k++)
                            {
                              f[i] += coeff * hist->f[i * 300 / nBin + k];
                              df[i] +=
                                coeff * coeff * hist->ff[i * 300 / nBin + k];
                            }
                          df[i] = sqrt (df[i] - f[i] * f[i] / hist->nPoints);
                        }

                      if (nin_ == 2)
                        strcpy (yname, "Diff. cross section [pb");
                      else
                        strcpy (yname, "Diff. width [GeV");
                      xName (hist->key, hist->plist, xname, units);
                      if (hist->restfr)
                        {
                          ii = 0;
                          while (((hist->restfr[ii]) != 0) && (ii < 10))
                            {
                              restnumbers[ii] = hist->restfr[ii] + '0';
                              ii++;
                            }
                          if (ii != 0)
                            strcat (xname, ", rest frame ");
                          while (ii < 10)
                            {
                              restnumbers[ii] = ' ';
                              ii++;
                            }
                          trim (restnumbers);
                          strcat (xname, restnumbers);
                        }
                      if (units[0])
                        {
                          strcat (yname, "/");
                          strcat (yname, units);
                        }
                      strcat (yname, "]");

                      plot_0 (hist->hMin, hist->hMax, nBin, f, df,
                              proces_1.proces, xname, yname);
                    }
              }
          }
        }
    }
}

void
editHist (void)
{
  do
    edittable (1, 4, &histTab, 1, "n_distrib", 0);
  while (correctHistList ());
}