#include "TROOT.h"
#include "gpData.h"
#include "TH1.h"
#include "TFile.h"
#include <iostream>
#include <string>

int main(int argc,char **argv) {

  string inputpath;
  Int_t beamenergy;
  string inputname;

  if(argc == 3) {
    inputpath = string(argv[1]);
    beamenergy = atoi(argv[2]);
  }
  else {
    cout << "usage: GPAnal_fz inputpath beam_energy" << endl;
    return 1;
  }

  Int_t start = inputpath.find_last_of("/");
  string outputname = inputpath.substr(start+1,inputpath.length())+ ".root";

  cout << "\nbeam energy: " << beamenergy << "\n";
  cout << "path is " << inputpath << "\n\n";
  cout << "root filename is " << outputname << "\n";

  //get all data files needed

  gpData d((inputpath+"/pairs.dat").c_str(),"pair");
  //cout << "got pairs \n";
  d.loadFile((inputpath+"/hadrons.dat").c_str(),"hadron");
  d.loadFile((inputpath+"/beam1.dat").c_str(),"beam1");
  //cout << "got beam1 \n";
  d.loadFile((inputpath+"/beam2.dat").c_str(),"beam2");
  //cout << "got beam2 \n";
  d.loadFile((inputpath+"/photon.dat").c_str(),"photon");
  //cout << "got photon \n";
  d.loadFile((inputpath+"/lumi.dat").c_str(),"lumi");

  cout << "finished loading files \n";

  TFile *output = new TFile(outputname.c_str(),"RECREATE");
  cout << "created root object \n";

  //make histogram objects
  TH1D *pae = new TH1D("pae","pair creation energies",1000,0,1);

  TH1D *b1e = new TH1D("b1e","energy of electrons in beam 1",1000,0,1);
  TH1D *b2e = new TH1D("b2e","energy of positrons in beam 2",1000,0,1);

  TH1D *b1thx = new TH1D("b1thx","theta-x of electrons in beam 1",1000,-5E-4,5E-4);
  TH1D *b2thx = new TH1D("b2thx","theta-x of positrons in beam 2",1000,-5E-4,5E-4);

  TH1D *b1thy = new TH1D("b1thy","theta-y of electrons in beam 1",1000,-5E-4,5E-4);
  TH1D *b2thy = new TH1D("b2thy","theta-y of positrons in beam 2",1000,-5E-4,5E-4);

  TH1D *b1x = new TH1D("b1x","x of electrons in beam 1",1000,-1000 ,1000 );
  TH1D *b2x = new TH1D("b2x","x of positrons in beam 2",1000,-1000 ,1000  );

  TH1D *b1y = new TH1D("b1y","y of electrons in beam 1",1000,-100 ,100  );
  TH1D *b2y = new TH1D("b2y","y of positrons in beam 2",1000,-100 ,100 );


  TH1D *lpe = new TH1D("lpe","energy of lumi electrons",1000,0,1);
  TH1D *lpp = new TH1D("lpp","energy of lumi positrons",1000,0,1);

  TH1D *lpx = new TH1D("lpx","x of lumi electrons",1000,-1000,1000);
  TH1D *lpy = new TH1D("lpy","y of lumi electrons",1000,-100 ,100);
  TH1D *lpz = new TH1D("lpz","z of lumi electrons",1000,-1000 ,1000);

  TH1D *pax = new TH1D("pax","x of pair e+/e-",1000,-1000 ,1000);
  TH1D *pay = new TH1D("pay","y of pair e+/e-",1000,-100 ,100);
  TH1D *paz = new TH1D("paz","z of pair e+/e-",1000,-1000 ,1000);

  TH1D *phe = new TH1D("phe","energy of photons",1000,0,1);


  //displaying some useful stats
  cout << "number of pairs: "<< d.numOfPcles("pair") << endl;
  cout << "number of hadrons: "<< d.numOfPcles("hadron") << endl;
  cout << "number of photons: "<< d.numOfPcles("photon") << endl;
  cout << "number of luminosity data points: "<< d.numOfPcles("lumi") << endl;
  cout << "number of macroparticles in beam 1: "<< d.numOfPcles("beam1") << endl;
  cout << "number of macroparticles in beam 2: "<< d.numOfPcles("beam2") << endl;


  //fill the histograms

  for(Int_t i=0;i<d.numOfPcles("pair");i++) {
    pair_part pp;
    d.getPcle(i,pp);

    pae->Fill(pp.e/beamenergy);
    pax->Fill(pp.x);
    pay->Fill(pp.y);
    paz->Fill(pp.z);

  }

  for(Int_t i=0;i<d.numOfPcles("beam1");i++) {
    beam1_part bp;
    d.getPcle(i,bp);

    b1thx->Fill(bp.tx);
    b1thy->Fill(bp.ty);
    b1y->Fill(bp.y);
    b1x->Fill(bp.x);
    b1e->Fill(bp.e/beamenergy);

  }

  for(Int_t i=0;i<d.numOfPcles("beam2");i++) {
    beam2_part bp;
    d.getPcle(i,bp);

    b2thx->Fill(bp.tx);
    b2thy->Fill(bp.ty);
    b2y->Fill(bp.y);
    b2x->Fill(bp.x);
    b2e->Fill(bp.e/beamenergy);

  }



  for(Int_t i=0;i<d.numOfPcles("lumi");i++) {
    lumi_part lp;
    d.getPcle(i,lp);

    lpe->Fill(lp.eele/beamenergy);
    lpp->Fill(lp.epos/beamenergy);
    lpx->Fill(lp.x);
    lpy->Fill(lp.y);
    lpz->Fill(lp.z);

  }

  for(Int_t i=0;i<d.numOfPcles("photon");i++) {
    photon_part p;
    d.getPcle(i,p);
    phe->Fill(p.e/beamenergy);
  }


  output->Write();
  output->Close();

  return 0;
}