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

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

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Program to plot PDF of Cherenkov light from muon using interpolation tables.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JAbstractHistogram<double> JHistogram_t;

  vector<string> inputFile;
  string         outputFile;
  double         E;
  double         R;
  JAngle3D       dir;
  double         TTS_ns;
  JHistogram_t   histogram;
  int            debug;

  try { 

    JParser<> zap("Program to plot PDF of Cherenkov light from muon using interpolation tables.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)                              = "";
    zap['E'] = make_field(E,         "muon energy [GeV]")          = 1.0;
    zap['R'] = make_field(R,         "distance of approach [m]");
    zap['D'] = make_field(dir,       "(theta, phi) of PMT [rad]");
    zap['T'] = make_field(TTS_ns,    "PMT time smearing [ns]")     = 0.0;   // [ns]
    zap['H'] = make_field(histogram, "histogram binning")          = JHistogram_t();
    zap['d'] = make_field(debug)                                   = 0;

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


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

  const int N = inputFile.size();

  int    type[N];
  JPDF_t pdf [N];

  try {

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

      NOTICE("loading input from file " << inputFile[i] << "... " << flush);
      
      type[i] = getPDFType(inputFile[i]);

      pdf [i].load(inputFile[i].c_str());

      pdf [i].setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
      
      if (TTS_ns > 0.0) {
	pdf[i].blur(TTS_ns);
      }
      
      NOTICE("OK" << endl);
    }
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }


  if (outputFile == "") {

    for (double dt; cin >> dt; ) {
      
      for (int i = 0; i != N; ++i) {

	JFunction1D_t::result_type y = pdf[i](R, dir.getTheta(), dir.getPhi(), dt);

	if        (is_bremsstrahlung(type[i])) {
	  y *= E;
	} else if (is_deltarays(type[i])) {
	  y *= getDeltaRaysFromMuon(E);
	}
      
	cout << setw(2)         << type[i]        << ' '
	     << SCIENTIFIC(7,1) << E              << ' '
	     << FIXED(5,1)      << R              << ' '
	     << FIXED(5,2)      << dir.getTheta() << ' '
	     << FIXED(5,2)      << dir.getPhi()   << ' '
	     << FIXED(5,1)      << dt             << ' '
	     << SCIENTIFIC(9,3) << get_value(y)   << endl;
      }
    }

    return 0;
  }


  TFile out(outputFile.c_str(), "recreate");
    
  int function = 0;
    
  if (inputFile.size() == 1 && getPDFType(inputFile[0]) == DIRECT_LIGHT_FROM_MUON) {
    function = 1;
  }
    
  //const double z0   = -100.0;                                 // [m]
  //const double t0   = (-z0 +  R * getTanThetaC()) / C;        // time [ns]
  const double t0   =  0.0;                                   // time [ns]

  if (!histogram.is_valid()) {

    if (function == 1) {

      histogram = JHistogram_t(t0 - 20.0, t0 + 50.0);

      histogram.setBinWidth(0.1);

    } else {

      histogram = JHistogram_t(t0 - 20.0, t0 + 500.0);

      histogram.setBinWidth(0.5);
    }
  }
  
  TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
  TH1D h1("h1", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
  TH1D h2("h2", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
    
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
      
    const double dt = h0.GetBinCenter(i) - t0;

    JFunction1D_t::result_type Y = JMATH::zero;

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

      JFunction1D_t::result_type y = pdf[j](R, dir.getTheta(), dir.getPhi(), dt);

      if        (is_bremsstrahlung(type[j])) {
	y *= E;
      } else if (is_deltarays(type[j])) {
	y *= getDeltaRaysFromMuon(E);
      }

      Y += y;
    }      

    h0.SetBinContent(i, get_value     (Y));
    h1.SetBinContent(i, get_derivative(Y));
    h2.SetBinContent(i, get_integral  (Y));
  }

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