#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 "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JQuantiles.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 * Demonstration program to plot RMS of arrival time of first hit as a function of
 * the minimal distance of approach of a muon to the PMT.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

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

  try { 

    JParser<> zap("Demonstration program to plot RMS of arrival time of first hit as a function of the minimal distance of approach of a muon to the PMT.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)       = "demo.root";
    zap['E'] = make_field(E_GeV,     "muon energy [GeV]");
    zap['R'] = make_field(R_Hz,      "background rate [Hz]");
    zap['D'] = make_field(dir,       "(theta, phi) of PMT [rad]");
    zap['d'] = make_field(debug)            = 0;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  typedef JSplineFunction1S_t                                     JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
		   JPolint1FunctionalGridMap,
		   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(zero));

  const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
			       SCATTERED_LIGHT_FROM_MUON,
			       DIRECT_LIGHT_FROM_EMSHOWERS,
			       SCATTERED_LIGHT_FROM_EMSHOWERS };
    
  const int  N = sizeof(pdf_t) / sizeof(pdf_t[0]);

  JPDF_t pdf[N];

  for (int i = 0; i != N; ++i) {

    try {

      const string file_name = getFilename(inputFile, pdf_t[i]);
      
      NOTICE("loading input from file " << file_name << "... " << flush);
      
      pdf[i].load(file_name.c_str());
      
      NOTICE("OK" << endl);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
    
    pdf[i].setExceptionHandler(supervisor);
  }


  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] = pdf[k](R, dir.getTheta(), dir.getPhi(), 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].f         +
	p[1].f         +
	p[2].f * E_GeV +
	p[3].f * 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);

    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 > notice_t) {
    timer.print(cout, w);
  }

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