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

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

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JPhysics/JPDF.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JTools/JElement.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JQuantiles.hh"


/**
 * Scaling of absorption and scattering length.
 */
double  absorptionLengthFactor;
double  scatteringLengthFactor;

inline double getAbsorptionLength(const double lambda)
{
  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
}


inline double getScatteringLength(const double lambda)
{
  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
}


/**
 * Use attenutation.
 */
bool    useAttenuation;


/**
 * Get attenuation length.
 *
 * \param  lambda   wavelength of photon [nm]
 * \return          attenuation length [m]
 */
inline double getAttenuationLength(const double lambda)
{
  static const double f = 0.17;

  const double l_abs = getAbsorptionLength(lambda);
  const double ls    = getScatteringLength(lambda) / f;

  double y = 1.0 / l_abs;

  if (useAttenuation) 
    y +=  1.0 / ls;

  return 1.0 / y;
}


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 draw PDF.
 */
int main(int argc, char **argv)
{
  using namespace std;

  string      outputFile;
  int         numberOfPoints;
  double      epsilon;
  double      E;
  double      R;
  orientation dir;
  vector<int> function;
  int         debug;

  try { 

    JParser<> zap;

    zap['o'] = make_field(outputFile)             = "pdf.root";
    zap['n'] = make_field(numberOfPoints)         = 25;
    zap['e'] = make_field(epsilon)                = 1.0e-10;
    zap['A'] = make_field(absorptionLengthFactor) = 1.0;
    zap['S'] = make_field(scatteringLengthFactor) = 1.0;
    zap['E'] = make_field(E)                      = 1.0;
    zap['R'] = make_field(R);
    zap['D'] = make_field(dir);
    zap['F'] = make_field(function);
    zap['U'] = make_field(useAttenuation);
    zap['d'] = make_field(debug)                  = 0;

    zap['F'] = JPARSER::not_initialised();

    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;

  // set global parameters

  const double P_atm = 300.0;       // ambient pressure [Atm]
  const double wmin  = 320.0;       // minimal wavelength [nm]
  const double wmax  = 610.0;       // maximal wavelength [nm]

  const JPDF_C
    pdf(NAMESPACE::getPhotocathodeArea(),
	NAMESPACE::getQE,
	NAMESPACE::getAngularAcceptance,
	getAttenuationLength,
	getScatteringLength,
	NAMESPACE::p00075,
	P_atm,
	wmin,
	wmax,
	numberOfPoints,
	epsilon);
  
 
  TFile out(outputFile.c_str(), "recreate");

  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 -   20.0;
  double xmax = t0 +  500.0;
  int    nx   = (int) ((xmax - xmin) / 0.5);

  if (function.size() == 1 && function[0] == 1) {
    xmin = t0 -   20.0;
    xmax = t0 +   50.0;
    nx   = (int) ((xmax - xmin) / 0.1);
  }
  
  TH1D h0("h0", NULL, nx, xmin, xmax);

  JSplineFunction1S_t f1;
 
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {

    const double dt = h0.GetBinCenter(i) - t0;

    double value = 0.0;

    for (vector<int>::const_iterator j = function.begin(); j != function.end(); ++j) {

      switch (*j) {

      case 1:
	value += pdf.getDirectLightFromMuon        (R, theta, phi, dt);
	break;

      case 2:
	value += pdf.getScatteredLightFromMuon     (R, theta, phi, dt);
	break;

      case 3:
	value += pdf.getDirectLightFromEMshowers   (R, theta, phi, dt) * E;
	break;
	
      case 4:
	value += pdf.getScatteredLightFromEMshowers(R, theta, phi, dt) * E;
	break;

      default:
	ERROR("Illegal function " << *j << endl);
	break;
      }
    }

    h0.SetBinContent(i, value);

    f1[dt] = value;
  }

  f1.compile();

  const double T = 5;  // [ns]

  JQuantiles quantiles(f1);

  const double t1 = quantiles.getX();
  const double y  = f1(t1 + T).v - f1(t1 - T).v;

  DEBUG("E     " << E                       << endl);
  DEBUG("R     " << R                       << endl);
  DEBUG("theta " << theta                   << endl);
  DEBUG("phi   " << phi                     << endl);
  DEBUG("int   " << quantiles.getIntegral() << endl);
  DEBUG("t1    " << t1                      << endl);
  DEBUG("max   " << quantiles.getY()        << endl);
  DEBUG("FWHM  " << quantiles.getFWHM()     << endl);
  DEBUG("int[] " << y                       << endl);  

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