#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <assert.h>

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

#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JQuantiles.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 determine PDF of first hit.
 */
int main(int argc, char **argv)
{
  using namespace std;

  string      inputFile;
  string      outputFile;  
  double      E_GeV;
  double      R_Hz;
  orientation dir;
  int         debug;

  try { 

    JParser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)       = "demo.root";
    zap['E'] = make_field(E_GeV);
    zap['R'] = make_field(R_Hz);
    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);
  }
  

  using namespace JTOOLS;
  using namespace JPHYSICS;
  using namespace JEEP;

  typedef JSplineFunction1S_t                                     JFunction1D_t;
  typedef
    JMapList<JPolint1FunctionalMap,
    JMapList<JPolint1FunctionalGridMap,
    JMapList<JPolint1FunctionalGridMap> > >                       JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  JPDF_t f[4];


  JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(0.0));
    
  for (int k = 0; k != 4; ++k) {

    string file_name = inputFile;

    string::size_type j = file_name.find('%');

    if (j != string::npos) {
    
      file_name[j] = '1' + k;
      
      try {

	NOTICE("loading input from file " << file_name << "... " << flush);
      
	f[k].load(file_name.c_str());

	NOTICE("OK" << endl);
      }
      catch(const JException& error) {
	FATAL(error.what() << endl);
      }
      
      for (JPDF_t::super_iterator i = f[k].super_begin(); i != f[k].super_end(); ++i)
	(*i).getValue().setExceptionHandler(supervisor);
    }
  }


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


  vector<double>           X;        // transverse distance [m]

  if (X.empty()) {
    X.push_back(  5.0);
    X.push_back( 15.0);
    X.push_back( 25.0);
    X.push_back( 35.0);
    X.push_back( 45.0);
    X.push_back( 55.0);
    X.push_back( 65.0);
    X.push_back( 75.0);
    X.push_back( 85.0);
    X.push_back( 95.0);
  }
  

  TFile out(outputFile.c_str(), "recreate");

  vector<TH1D*> buffer;

  for (vector<double>::const_iterator i = X.begin(); i != X.end(); ++i) {

    ostringstream os;

    os << "t[" << *i << "]";
    
    buffer.push_back(new TH1D(os.str().c_str(), NULL, 5000, -250.0, +250.0));
  }

  vector<double> x;

  {
    vector<double>::const_iterator j = X.begin();
    vector<double>::const_iterator i = j++;

    x.push_back(*i - 0.5*(*j-*i));

    for (; j != X.end(); ++i, ++j) 
      x.push_back(0.5*(*i+*j));

    --i;
    --j;

    x.push_back(*j + 0.5*(*j-*i));
  }

  TH1D rms("rms", NULL, X.size(), &x[0]);


  JFunction1D_t::result_type p[4];

  JTimer timer;

  timer.reset();

  const double Tmin = -250.0;     // [ns]
  const double Tmax = +250.0;     // [ns]

  int n = 0;

  for (size_t i = 0; i != X.size(); ++i) {

    TH1D*        h = buffer[i];
    const double R = X[i];

    JSplineFunction1S_t zsp;

    for (int j = 1; j <= h->GetNbinsX(); ++j) {

      const double t1 = h->GetBinCenter(j);

      timer.start();

      ++n;

      for (int k = 0; k != 4; ++k)
	p[k] = f[k](R, theta, phi, t1);

      // integral [Tmin,t1]
      const double v = 
	p[0].v         +
	p[1].v         +
	p[2].v * E_GeV +
	p[3].v * E_GeV +
	R_Hz * 1e-9 * (t1   - Tmin);

      // integral [Tmin,Tmax]
      const double V = 
	p[0].V         +
	p[1].V         +
	p[2].V * E_GeV +
	p[3].V * E_GeV +
	R_Hz * 1e-9 * (Tmax - Tmin);

      // function value
      const double y = 
	p[0].y         +
	p[1].y         +
	p[2].y * E_GeV +
	p[3].y * E_GeV + 
	R_Hz * 1e-9;

      const double W = exp(-v) * y / (1.0 - exp(-V));

      timer.stop();

      zsp[t1] = W;

      h->SetBinContent(j,W);
    } 

    zsp.compile();

    JQuantiles result(zsp, 1.0e-10);

    const double sig = result.getFWHM() * 0.5 / sqrt(2.0*log(2.0));

    rms.Fill(R, sig);

    NOTICE("integral " << R << ' ' << result.getIntegral() << endl);
  }

  const float w = 1.0 / (float) n;

  if (debug > 1)
    timer.print(cout, w);
  
  out.Write();
  out.Close();
}