#include "HepTuple/HepHBookFileManager.h"
#include "ZMutility/iostream"
#ifndef ISOCXX__ISOCXX
  #include <STLUtility/fwd_vector.h>
#endif
#include <vector>
#include <string>
#include "HepTuple/HepNtuple.h"
#include "HepTuple/HepHist1D.h"
USING( namespace std )
#define MAXTRACKS 200

int main(){

/******* Initialise hbook files *******/
  //Open file and connect to input ntuple.
  HepFileManager *outManager;
  outManager = new HepHBookFileManager("new.hbk");

  HepFileManager *inManager;
  inManager = new HepHBookFileManager("sample.hbk");

  HepNtuple &inTuple = inManager->retrieveNtuple("Test N-tuple");
  inTuple.setColumnWise();
  inTuple.setDiskResident();

  //Open new histogram for output and define columns
  HepNtuple &outTuple = outManager->ntuple("Smeared",20);
  outTuple.setColumnWise();
  outTuple.setDiskResident();
  outTuple.columnArray("Track::IP",1,(Float8)0.);
  outTuple.dimension("Track::IP",MAXTRACKS);
  outTuple.columnArray("Track::Phi0",1,(Float8)0.);
  outTuple.dimension("Track::Phi0",MAXTRACKS);
  outTuple.columnArray("Track::Curv",1,(Float8)0.);
  outTuple.dimension("Track::Curv",MAXTRACKS);
  outTuple.clearData();

  outManager->mkdir("Histograms");
  outManager->cd("Histograms");

  HepHist1D &outHist = outManager->hist1D("ntracks",100,-100,100.0);

 //  Int4 nevts =  inTuple.nrows();
  Int4 nevts = 10;
  std::cout << "File contains " << nevts << " events.\n\n";

  for(int n=0; n<nevts;) {    // Loop over events
    ++n;

/******* Read data from ntuple *******/
    Int4 ntracks;
    Int4 idpart[MAXTRACKS];
    Float4 mom[MAXTRACKS][5];
    Float4 vtx[MAXTRACKS][5];
    Int4 p[MAXTRACKS];
    Int4 lp[MAXTRACKS];
    inTuple.destinationAt("OURVARS1::npart",&ntracks);
    inTuple.destinationAt("OURVARS1::idpart",(Int4*)&(idpart[0]));
    inTuple.destinationAt("OURVARS1::Pmom",&(mom[0][0]));
    inTuple.destinationAt("OURVARS1::Vtx",&(vtx[0][0]));
    inTuple.destinationAt("OURVARS1::IDparent",&(p[0]));
    inTuple.destinationAt("OURVARS1::linePar",&(lp[0]));
    inTuple.readBlock(n,"OURVARS1");
    ntracks = MAXTRACKS;

    std::cout << "Event " << n << " has " << ntracks << " tracks.\n";

    Float8 ip[MAXTRACKS],phi[MAXTRACKS],curv[MAXTRACKS];
    // Try loading some variables for later use (depends on data structure)**
    int i,j;
    for (i=0; i<ntracks; ++i){
      ip[i] = mom[i][1] = 0.0;
      phi[i] = mom[i][2]= 0.0;
      curv[i] = mom[i][3]= 0.0;
    }
    for (;i<MAXTRACKS;++i){
      ip[i] = 0.0;
      phi[i] = 0.0;
      curv[i] = 0.0;
    }

/******* Write data to ntuple *******/
    for (j=0; j<ntracks; ++j){
      for (i=0; i<5; ++i){
	outHist.accumulate(vtx[i][j]);
      }
    }

    outTuple.capture("Track::IP",ip);
    outTuple.capture("Track::Phi0",phi);
    outTuple.capture("Track::Curv",curv);

    outTuple.storeCapturedData();
    outTuple.clearData();


  } // *Event Loop*

/******* Close hbook files etc *******/
  //inManager->write();
  outManager->write();

  delete outManager;
  delete inManager;  // implicit HROUT

  return 0;
}