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

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

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JTransformer.hh"
#include "JPhysics/JPDFTable.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; }


/**
 * Auxiliary program to plot PDF.
 */
int main(int argc, char **argv)
{
  using namespace std;

  string      inputFile;
  string      outputFile;
  double      D;
  double      cd;
  orientation dir;
  int         debug;

  try { 

    JParser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)       = "";
    zap['R'] = make_field(D);
    zap['c'] = make_field(cd);
    zap['D'] = make_field(dir);
    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 JMultipleMap<4, JPolint1FunctionalMap>::type_list       JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  JPDF_t pdf;


  try {

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

    pdf.load(inputFile.c_str());

    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
    
  JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(0.0));
    
  for (JPDF_t::super_iterator i = pdf.super_begin(); i != pdf.super_end(); ++i)
    (*i).getValue().setExceptionHandler(supervisor);


  if (outputFile == "") {

    for ( ; ; ) {
    
      double  dt;
      
      cout << "> " << flush;
      cin  >> dt;
      
      const double value = get_value(pdf(D, cd, theta, phi, dt));
      
      cout << endl;
      cout << value << endl;
    }

  } else {


    TFile out(outputFile.c_str(), "recreate");
    
    int function = 13;
    
    if (inputFile.find("13") == string::npos) 
      function = 0;
    
    const double t0   =  D * getIndexOfRefraction() / C;        // time [ns]
    //const double t0   =  0.0;                                   // time [ns]

    double xmin = t0 -   10.0;
    double xmax = t0 + 1000.0;
    int    nx   = (int) ((xmax - xmin) / 0.5);
    
    if (function == 13) {
      xmin = t0 -  10.0;
      xmax = t0 +  20.0;
      nx   = (int) ((xmax - xmin) / 0.1);
    }
    
    TH1D h0("h0", NULL, nx, xmin, xmax);
    
    for (int i = 1; i <= h0.GetNbinsX(); ++i) {
      
      const double dt = h0.GetBinCenter(i) - t0;

      const double value = get_value(pdf(D, cd, theta, phi, dt));
      
      h0.SetBinContent(i, value);
    }

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