#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <map>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TRandom3.h"

#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JTransformer.hh"
#include "JPhysics/JCDFTable.hh"


typedef std::pair<double,double> orientation;

std::istream& operator>>(std::istream& in, orientation& x) { return in >> x.first >> x.second; }
std::ostream& operator<<(std::ostream& out, const orientation& x) { return out << x.first << ' ' << x.second; }


/**
 * Main program to verify hit generation
 */
int main(int argc, char **argv)
{
  using namespace std;

  string      inputFile;
  string      outputFile;
  double      E;
  double      R;
  orientation dir;
  int         numberOfEvents;
  int         debug;

  try { 

    JParser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)     = "";
    zap['E'] = make_field(E)              = 1.0;         // [GeV]
    zap['R'] = make_field(R);                            // [m]
    zap['D'] = make_field(dir);                          // [rad]
    zap['N'] = make_field(numberOfEvents) = 0;
    zap['d'] = make_field(debug)          = 0;

    if (zap.read(argc, argv) != 0)
      return 1;
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  const double theta = dir.first;
  const double phi   = dir.second;
  

  using namespace JPHYSICS;
  using namespace JTOOLS;


  typedef JSplineFunction1D_t                                     JFunction1D_t;
  //typedef JPolint1Function1D_t                                    JFunction1D_t;
  typedef
    JMapList<JPolint1FunctionalMap,
    JMapList<JPolint1FunctionalGridMap,
    JMapList<JPolint1FunctionalGridMap> > >                       JMapList_t;
  typedef JCDFTable<JFunction1D_t, JMapList_t>                    JCDF_t;

  JCDF_t cdf;


  try {

    NOTICE("loading input from file " << inputFile << "... " << flush);

    cdf.load(inputFile.c_str());

    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }


  if (outputFile == "") {

    for ( ; ; ) {

      double  x;

      cout << "> " << flush;
      cin  >> x;

      try {
        
        const double npe = cdf.getNPE (R, theta, phi) * E;        
        const double t   = cdf.getTime(R, theta, phi, x);

        cout << " --> " << t << ' ' << npe;
      }
      catch(const exception& error) {
        ERROR(error.what() << endl);
      }

      cout << endl;
    }

  } else {
  
   
    TFile out(outputFile.c_str(), "recreate");
  
    int function = 1;

    if (inputFile.find('1') == string::npos)
      function = 0;

    const double z0   = -100.0;                                 // [m]
    //const double t0   = (-z0 +  R * getTanThetaC()) / C;        // time [ns]
    const double t0   =  0.0;                                   // time [ns]

    double xmin = t0 -   10.0;
    double xmax = t0 +  500.0;
    int    nx   = (int) ((xmax - xmin) / 0.5);
  
    if (function == 1) {
      xmin = t0 -   20.0;
      xmax = t0 +   20.0;
      nx   = (int) ((xmax - xmin) / 0.1);
    }

    TH1D h0("h0", NULL, nx, xmin, xmax);
    TH1D h1("h1", NULL, 100000, 0.0, +1.0);

    h0.Sumw2();
  
    for (int j = 1; j <= h1.GetNbinsX(); ++j) {
    
      try {
      
	const double x = h1.GetBinCenter(j);
	const double t = cdf.getTime(R, theta, phi, x);
      
	h1.SetBinContent(j, t);
      }
      catch(const exception& error) {
	ERROR(error.what() << endl);
      }
    }
    
    if (numberOfEvents > 0) {
    
      using namespace JEEP;
  
      JTimer timer;
    
      timer.reset();
    
      for (int i = 0; i != numberOfEvents; ++i) {
      
	if (i%1000== 0) { NOTICE(i << '\r'); }
      
	timer.start();
      
	try {
	
	  const double npe = cdf.getNPE(R, theta, phi) * E;
	  
	  for (int j = gRandom->Poisson(npe); j != 0; --j) {
	  
	    const double x = gRandom->Rndm();
	    const double t = cdf.getTime(R, theta, phi, x);
	  
	    h0.Fill(t + t0);
	  }
	}
	catch(const exception& error) {
	  NOTICE(error.what() << endl);
	}
      
	timer.stop();
      }
      NOTICE(endl);
  
      const double w = 1.0 / (double) numberOfEvents;
    
      if (debug > 1)
	timer.print(cout, w);
    
      // normalise histogram (dP/dt)
    
      for (int j = 1; j <= h0.GetNbinsX(); ++j) {
      
	const double y  = h0.GetBinContent(j);
	const double z  = h0.GetBinError  (j);
	const double dx = h0.GetBinWidth  (j);
      
	h0.SetBinContent(j, y / (dx*numberOfEvents*E));
	h0.SetBinError  (j, z / (dx*numberOfEvents*E));
      }
    }

    out.Write();
    out.Close();
  }
}