#include <math.h>
#include <iostream.h>
#include "cfortran/cfortran.h"
#include "cfortran/hbook.h"
#include "ThreeVector.h"
#include "FourVector.h"
#include "anl_utils.h"

#define PAWC_SIZE 10000000
float pawc_[PAWC_SIZE];

void main()
{

// set up HBOOK, open ntuple file

  HLIMIT(PAWC_SIZE);
  int lun = 20;
  int istat = 0;
  int lrec = 1024;
  char chopt[2] = " ";
  char infile[81] = "evtanl.hbook";
  cout << "Enter ntuple file name" << endl;
  cin.getline (infile,80);
  HROPEN (lun, "ntuple", infile, chopt, lrec, istat);
  if (istat != 0){
    cout << "HROPEN error, istat = " << istat << endl;
    exit(istat);
  }

// read in ntuple only (here, ntupleId = 1).

  int ntupleId = 1;
  HRIN (ntupleId, 9999, 0);

// Book analysis histograms starting with 1000

  HBOOK1 (1000, "px1", 100, -100., 100., 0.);
  HBOOK1 (1001, "py1", 100, -100., 100., 0.);
  HBOOK1 (1002, "pz1", 100, -100., 100., 0.);
  HBOOK1 (1003, "E1", 100, 0., 200., 0.);
  HBOOK1 (1004, "btag1", 100, 0., 1., 0.);

  HBOOK1 (1005, "px2", 100, -100., 100., 0.);
  HBOOK1 (1006, "py2", 100, -100., 100., 0.);
  HBOOK1 (1007, "pz2", 100, -100., 100., 0.);
  HBOOK1 (1008, "E2", 100, 0., 200., 0.);
  HBOOK1 (1009, "btag2", 100, 0., 1., 0.);

  HBOOK1 (1010, "px3", 100, -100., 100., 0.);
  HBOOK1 (1011, "py3", 100, -100., 100., 0.);
  HBOOK1 (1012, "pz3", 100, -100., 100., 0.);
  HBOOK1 (1013, "E3", 100, 0., 200., 0.);
  HBOOK1 (1014, "btag3", 100, 0., 1., 0.);

  HBOOK1 (1015, "px4", 100, -100., 100., 0.);
  HBOOK1 (1016, "py4", 100, -100., 100., 0.);
  HBOOK1 (1017, "pz4", 100, -100., 100., 0.);
  HBOOK1 (1018, "E4", 100, 0., 200., 0.);
  HBOOK1 (1019, "btag4", 100, 0., 1., 0.);

// get information on ntuple

  int ncol;
  char title[80];
  const int MAX_COL = 50;
  float rlow[MAX_COL];
  float rhigh[MAX_COL];
  char tag[MAX_COL][8];
  HGIVEN (ntupleId, title, ncol, tag, rlow, rhigh);
  cout << "number of columns = " << ncol << endl;
  cout << "event variables..." << endl;
  for (int i=0; i<ncol; i++){ cout << i << "   " << tag[i] << endl; }

//  The analysis part comes here.  N.B. loop over events starts at 1, not 0.

  float* event = new float[ncol];
  FourVector pjet[4];
  float btag[4];
  int numEvents, ierr; 
  HNOENT(ntupleId, numEvents);
  cout << "number of events = " << numEvents << endl;
  char chrout[60] = "evtanl";
  HGNPAR(ntupleId, chrout);
  for (int i=1; i<=numEvents; i++){
    HGNF(ntupleId, i, *event, ierr);
    for (int j=0; j<ncol; j++){
      HF1 (1000 + j, event[j], 1.);
    }

    for (int j=0; j<4; j++){
      pjet[j] = pjet_from_event(j, event);
      btag[j] = event[j*5 + 4];
    }

    // put your analysis code here





  }
  HREND ("ntuple");             // close ntuple file






// Delete ntuple and store analysis histograms

  HDELET (ntupleId);
  int outlun = 30;
  char outfile[81] = "evtanl.hbook";
  cout << "Enter output histogram file name" << endl;
  cin.getline (outfile,80);
  HROPEN (outlun, "histog", outfile, "N", lrec, istat);
  if (istat != 0){
    cout << "HROPEN error, istat = " << istat << endl;
    exit(istat);
  }
  int icycle = 0;
  HROUT (0, icycle, " ");        // store histograms
  HREND ("histog");              // close histogram file

}