#include <iostream.h>
#include "cfortran/cfortran.h"
#include "cfortran/hbook.h"
#include "Utils.h"
#include "GenEvent.h"
#include "SimEvent.h"
#include "random.h"
#include "toymc.h"

#define PAWC_SIZE 10000000
float pawc_[PAWC_SIZE];

void main()
{

  const float pi = 3.14159265;
  const float mz = 91.2;

// set up HBOOK 

  HLIMIT(PAWC_SIZE);
  int lun = 20;
  int istat = 0;
  int lrec = 1024;
  char outfile[81] = "toymc.his";
  cout << "Enter histogram file name" << endl;
  cin.getline (outfile,80);
  HROPEN (lun, "histog", outfile, "N", lrec, istat);
  if (istat != 0){
    cout << "HROPEN error, istat = " << istat << endl;
    exit(istat);
  }

  const int Nrows = 20; 
  float evt[Nrows];
  char chtags[20][6] = {"px1", "py1", "pz1", "E1", "btag1",
                        "px2", "py2", "pz2", "E2", "btag2",
                        "px3", "py3", "pz3", "E3", "btag3",
                        "px4", "py4", "pz4", "E4", "btag4"};
                        
  HBOOKN (1, "event data", 20, " ", 10000, chtags);

  HBOOK1 (101, "phi", 100, 0., 2.*pi, 0.);
  HBOOK1 (102, "ct", 100, -1., 1., 0.);
  HBOOK1 (103, "phi_star", 100, 0., 2.*pi, 0.);
  HBOOK1 (104, "cts", 100, -1., 1., 0.);
  HBOOK1 (106, "mzstar", 100, 80., 100., 0.);

// set up parameters of simulation 

  int num_events;
  float Ecm, mh;
  cout << "Enter number of events" << endl;
  cin >> num_events;
  cout << "Enter Ecm" << endl;
  cin >> Ecm;

  int evtType;
  cout << "e+e- -> HZ ..............1" << endl;
  cout << "e+e- -> ZZ...............2" << endl;
  cout << "Enter event type" << endl;
  cin >> evtType;

  // Find fmax values 

  int seed = 12345;
  float fmax_u, fmax_d, fmax_ZZ;
  if ( evtType == 1 ) {
    cout << "Enter Higgs mass" << endl;
    cin >> mh;
    fmax_u = get_fmax_HZ (&seed, 100000, 1, Ecm, mh);
    cout << "fmax_u  = " << fmax_u << endl;
    fmax_d = get_fmax_HZ (&seed, 100000, 2, Ecm, mh);
    cout << "fmax_d  = " << fmax_d << endl;
    fmax_u  = fmax_u * 1.1;
    fmax_d  = fmax_d * 1.1;
  }
  else if ( evtType == 2) {
    fmax_ZZ = get_fmax_ZZ(&seed, 100000, Ecm, mz);
    cout << "fmax_ZZ = " << fmax_ZZ << endl;
    fmax_ZZ = fmax_ZZ * 1.1;
  }

  //  generate events

  for(int i = 1; i <= num_events; i++){

    GenEvent genEvt;
    if ( evtType == 1 ) {
      genEvt = generate_HZ (&seed, Ecm, mh, fmax_u, fmax_d); 
    }
    else if ( evtType == 2 ) {
      genEvt = generate_ZZ (&seed, Ecm, fmax_ZZ); 
    }

    SimEvent simEvt = simulate_det (&seed, genEvt);

    // pack event into array and fill ntuple

    evt[0] = simEvt.P1().x();
    evt[1] = simEvt.P1().y();
    evt[2] = simEvt.P1().z();
    evt[3] = simEvt.P1().t();
    evt[4] = simEvt.btag1();

    evt[5] = simEvt.P2().x();
    evt[6] = simEvt.P2().y();
    evt[7] = simEvt.P2().z();
    evt[8] = simEvt.P2().t();
    evt[9] = simEvt.btag2();

    evt[10] = simEvt.P3().x();
    evt[11] = simEvt.P3().y();
    evt[12] = simEvt.P3().z();
    evt[13] = simEvt.P3().t();
    evt[14] = simEvt.btag3();

    evt[15] = simEvt.P4().x();
    evt[16] = simEvt.P4().y();
    evt[17] = simEvt.P4().z();
    evt[18] = simEvt.P4().t();
    evt[19] = simEvt.btag4();

    HFN (1, evt);

  }

// Store histograms, ntuples and close

  int icycle=0;
  HROUT (0, icycle, " ");
  HREND ("histog");

}