//Author: Jordan Seneca
//Use: ./HistRecoTrks inpath* outfile

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <glob.h>
//#include "/sps/km3net/users/jseneca/docopt.cpp/docopt.h"

#include "util.hh"
#include "Vec.hh"
#include "Mat.hh"
#include "Trk.hh"
#include "Evt.hh"
#include "EventFile.hh"

#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TString.h"
#include "TTree.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "TFile.h"
#include "TCanvas.h"

#include "HistRecoTrks.hh"

static const char USAGE[] = 
R"(

HistRecoTrks.

    Usage:
        HistRecoTrks <indir> <outfile> <detfile>

)";


using namespace std;

int main(int argc, char *argv[] ){
  
  /*
  map<string, docopt::value> args = docopt::docopt(USAGE,
						   { argv + 1, argv + argc },
						   true,
						   "HistRecoTrks 0.1");

  for(auto const& arg: args){
    cout<< arg.first << arg.second << endl;
  }
  */

  cout << "File path to be searched: " << argv[1] << endl ;
  
  TFile outfile(argv[2], "RECREATE", "outfile");
  
  string DETECTOR_FILE( argv[3] );

  std::cout << "Still alive here?" << std::endl;
  const uint region_n = 3;

  //All compared against Energy for TH2D and Energy + NHits for TH3D  
  TH3D* h3dAng[region_n];
  TH3D* h4dAng[region_n];
  TH3D* h4dAngOld[region_n];

  TH3D* h3dErec[region_n];
  TH3D* h4dErec[region_n];
  TH3D* h4dErecOld[region_n];

  TH2D* h3dE[region_n];
  TH2D* h4dE[region_n];
  TH2D* hMcE[region_n];

  TH3D* hPos[region_n];

  for( uint iregion = 0; iregion < region_n; iregion++ ){

    vector< Double_t > angbin_vec = get_log_vector(1e-3, 180, 1000);
    for (Int_t i=0;i<=1000;i++) {
      
    }
    const vector< Double_t > Ebin_vec   = get_lin_vector(1e5,  1e8, 1000);
    const vector< Double_t > hitbin_vec = get_lin_vector(1e3,  1e5, 100);

    h3dAng[iregion] = new TH3D( (TString)"h3dAngENhits" + to_string(iregion), " ;delta ang[deg]; E[GeV]; hit_N", 1000, &angbin_vec[0], 1000, &Ebin_vec[0], 100, &hitbin_vec[0]);
    h4dAng[iregion] = new TH3D( (TString)"h4dAngENhits" + to_string(iregion), " ;delta ang[deg]; E[GeV]; hit_N", 1000, &angbin_vec[0], 1000, &Ebin_vec[0], 100, &hitbin_vec[0]);
    h4dAngOld[iregion] = new TH3D( (TString)"h4dAngOldENhits" + to_string(iregion), " ;delta ang[deg]; E[GeV]; hit_N", 1000, &angbin_vec[0], 1000, &Ebin_vec[0], 100, &hitbin_vec[0]);

    
    h3dErec[iregion] = new TH3D( (TString)"h3dErecENhits" + to_string(iregion), "E res comparison; delta E[deg]; E [GeV]; hit_N", 200, -2, 2, 1000, 1e5, 1e8, 100, 1e3, 1e5);
    h4dErec[iregion] = new TH3D( (TString)"h4dErecENhits" + to_string(iregion), "E res comparison; delta E[deg]; E [GeV]; hit_N", 200, -2, 2, 1000, 1e5, 1e8, 100, 1e3, 1e5);
    h4dErecOld[iregion] = new TH3D( (TString)"h4dErecOldENhits" + to_string(iregion), "E res comparison; delta E[deg]; E [GeV]; hit_N", 200, -2, 2, 1000, 1e5, 1e8, 100, 1e3, 1e5);


    h3dE[iregion] = new TH2D( (TString)"h3dENhits" + to_string(iregion), "3d E; E [GeV]; hit_N", 1000, 1e5, 1e8, 100, 1e3, 1e5);
    h4dE[iregion] = new TH2D( (TString)"h4dENhits" + to_string(iregion), "4d E; E [GeV]; hit_N", 1000, 1e5, 1e8, 100, 1e3, 1e5);
    hMcE[iregion] = new TH2D( (TString)"hMcENhits" + to_string(iregion), "mc E; E [GeV]; hit_N", 1000, 1e5, 1e8, 100, 1e3, 1e5);

    hPos[iregion] = new TH3D( (TString)"hPosition" + to_string(iregion), "pos mapping; pos x [m]; pos y [m]; pos z [m]", 300, -600, 600, 300, -600, 600, 400, 0, 800); //Don't make this histo too big or the file will reach the byte limit
    
  }

  TH2D* h3dAngDist = new TH2D("h3dAngDist", " ;delta ang[deg]; d_can [m]", 1000, 0, 180, 200, -300, 550);
  TH2D* h4dAngDist = new TH2D("h4dAngDist", " ;delta ang[deg]; d_can [m]", 1000, 0, 180, 200, -300, 550);

  TH2D* h3dErecDist = new TH2D("h3dErecDist", "E res comparison; delta E[deg]; d_can [m]", 200, -2, 2, 200, -300, 550);
  TH2D* h4dErecDist = new TH2D("h4dErecDist", "E res comparison; delta E[deg]; d_can [m]", 200, -2, 2, 200, -300, 550);

  //This against HitN
  TH2D* h3dStopwatch = new TH2D("h3dStopwatch", " ;reco_T [s]; hit_N", 1000, 0, 1e3, 1000, 0, 1e5);
  TH2D* h4dStopwatch = new TH2D("h4dStopwatch", " ;reco_T [s]; hit_N", 30,   0, 1e3, 1000, 0, 1e5);


  vector<string> filelist = {};

  glob_t globbuf;

  glob(argv[1], 0, NULL, &globbuf);
  for(uint i = 0; i < globbuf.gl_pathc; i++){
    filelist.push_back( globbuf.gl_pathv[i] );
    cout << "Added file: " << globbuf.gl_pathv[i] << endl;
  }

  cout << "Filelist length: " << filelist.size() << endl;

  Vec coord_origin;
  for(vector<string>::const_iterator filename = filelist.begin(); filename != filelist.end(); filename++){
    EventFile file( (*filename).c_str() );

    cout << "Found event file " <<  *filename << endl;
    
    Head header = file.header;
    Det detector( DETECTOR_FILE );
    DetDim detdim( detector );
    cout << "Height: " << detdim.height() << endl;
    cout << "Radius: " << detdim.radius() << endl;

    try{
      coord_origin = header.coord_origin();
    } catch(std::invalid_argument err){
      //cout << "Error, :" << err << endl;
      coord_origin = detector.get_cg();
    }

    
    uint n_events = 0;
    uint n_contained = 0;
    uint n_inside    = 0;
    uint n_outside   = 0;

    uint iregion = 0; //0 contained, 1 inside, 2 outside

    for( EventFile::iterator ievt = file.begin(); ievt != file.end(); ievt++ ){
      Evt evt = *ievt;
      n_events += 1;
      
      Trk mc_neut = evt.mc_trks[0];
      mc_neut.pos += coord_origin;

      const double dist_edge = get_dist_to_edge( mc_neut, detector, detdim );

      if( dist_edge < -100 ){
	//cout << "Event contained: " << dist_edge << endl;
	iregion = 0; n_contained++;
      } else if( dist_edge < 0 ){
	//cout << "Event inside: " << dist_edge << endl;
	iregion = 1; n_inside++;
      } else {
	//cout << "Event outside: " << dist_edge << endl;
	iregion = 2; n_outside++;
      };


      //cout << "evt.hits.size(): " << evt.hits.size() << endl;

      hMcE[iregion]->Fill( mc_neut.E, evt.hits.size() );
      hPos[iregion]->Fill( mc_neut.pos.x, mc_neut.pos.y, mc_neut.pos.z );
      
      for( vector<Trk>::iterator trk =evt.trks.begin(); trk != evt.trks.end(); trk ++ ){

	if(trk->rec_type == 103){
	  Reso myres( *trk, mc_neut);
	  h3dAng[iregion] ->Fill( myres.ang_res, mc_neut.E,  evt.hits.size() );
	  h3dErec[iregion]->Fill( myres.E_res  , mc_neut.E,  evt.hits.size() );
	  h3dE[iregion]   ->Fill(                   trk->E,  evt.hits.size() );
	  //h3dStopwatch    ->Fill( 0,                         evt.hits.size() );
	  h3dAngDist      ->Fill( myres.ang_res, dist_edge);
	  h3dErecDist     ->Fill( myres.E_res  , dist_edge);
	  continue;
	}
        
	if(trk->rec_type == 104){
	  Reso myres(*trk, mc_neut);
	  h4dAng[iregion] ->Fill( myres.ang_res, mc_neut.E,  evt.hits.size() );
	  h4dErec[iregion]->Fill( myres.E_res  , mc_neut.E,  evt.hits.size() );
	  h4dE[iregion]   ->Fill(                   trk->E,  evt.hits.size() );
	  h4dStopwatch    ->Fill( evt.getusr("newshowerfit_time"),      evt.hits.size() );
	  h4dAngDist      ->Fill( myres.ang_res, dist_edge);
	  h4dErecDist     ->Fill( myres.E_res  , dist_edge);
	  continue;
	}

	if(trk->rec_type == 114){
	  Reso myres(*trk, mc_neut);
	  h4dAngOld[iregion] ->Fill( myres.ang_res, mc_neut.E,                   evt.hits.size() );
	  h4dErecOld[iregion]->Fill( myres.E_res  , mc_neut.E,                   evt.hits.size() );
	  continue;
	}
      }
    }
    cout << "analysed "   << n_events    << " events" << endl;
    cout << "contained: " << n_contained << " events" << endl;
    cout << "inside   : " << n_inside    << " events" << endl;
    cout << "outside  : " << n_outside   << " events" << endl;
  }

  cout << "Writing file to " << argv[2] << "... " << endl;
  outfile.Write();
  cout << "Closing file... " << endl;
  outfile.Close();
  
  cout << "Hists done." << endl;
};