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

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiMapTransformer.hh"
#include "JTools/JQuadrature.hh"

#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JGeanx.hh"
#include "JPhysics/JGeanz.hh"

#include "JSystem/JSystemToolkit.hh"

#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

#include "JConstants.hh" 

/**
 * \file
 *
 * Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a shower.
 *
 * The PDFs are tabulated as a function of <tt>(E, D, cos(theta), theta, phi, t)</tt>, where:
 * - <tt>E</tt> is the logarithm of the energy;
 * - <tt>D</tt> is the distance between the vertex and the PMT;
 * - <tt>cos(theta)</tt> the cosine of the photon emission angle;
 * - <tt>(theta, phi)</tt> the orientation of the PMT; and
 * - <tt>t</tt> the arrival time of the light with respect to the Cherenkov hypothesis.
 *
 * The orientation of the PMT is defined in the coordinate system in which
 * the shower developes along the z-axis and the PMT is located in the x-z plane.
 * \author jseneca, mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string              fileDescriptor;
  string              outputFile;
  int                 numberOfPoints;
  size_t              elong_sample_count;
  size_t              cosine_angle_count; // Cosine angle is not as sharp when elongated
  double              epsilon;
  double              TTS;
  set<double>         Ds;                          // distance [m]
  set<double>         Es;                         // logarithm of energy (GeV)
  bool                option;
  int                 debug;
  
  try {

    JParser<> zap("Program to create interpolation tables of the PDF of the arrival time of the Cherenkov light from a shower.");

    zap['F'] = make_field(fileDescriptor,         "file name descriptor for PDF tables");
    zap['o'] = make_field(outputFile,             "file name output PDF");
    zap['n'] = make_field(numberOfPoints,         "points for integration")     = 25;
    zap['e'] = make_field(epsilon,                "precision for integration")  = 1.0e-10;
    zap['D'] = make_field(Ds,                     "distance [m]")               = JPARSER::initialised();
    zap['E'] = make_field(Es,                     "energy [log10(GeV)]")        = JPARSER::initialised();
    zap['C'] = make_field(cosine_angle_count,     "points for cosine angle")    = 60;
    zap['N'] = make_field(elong_sample_count,     "points for elongation")      = 20;
    zap['T'] = make_field(TTS,                    "sigma of time bluring [ns]") = 2.0;
    zap['O'] = make_field(option,                 "optional transformation");
    zap['d'] = make_field(debug)                  = 0;

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


  if (cosine_angle_count == 0) {
    FATAL("Invalid option -C <" << cosine_angle_count << ">" << endl);
  }

  typedef JPolint1Function1D_t                                    JFunction1D_t;
  
  typedef JMAPLIST<JPolint1FunctionalMap,
		   JPolint1FunctionalMap,
		   JPolint1FunctionalGridMap,
		   JPolint1FunctionalGridMap>::maplist            J4DMapList_t;

  typedef JMAPLIST<JPolint1FunctionalMap,
		   JPolint1FunctionalMap,
		   JPolint1FunctionalMap,
		   JPolint1FunctionalGridMap,
		   JPolint1FunctionalGridMap>::maplist            J5DMapList_t;

  typedef JPDFTable<JFunction1D_t, J4DMapList_t>                  JPDF4d_t;
  typedef JPDFTable<JFunction1D_t, J5DMapList_t>                  JPDF5d_t;

  typedef JPDFTransformer<4, JFunction1D_t::argument_type>        JFunction4DTransformer_t;
  typedef JPDFTransformer<5, JFunction1D_t::argument_type>        JFunction5DTransformer_t;
  typedef JArray<5, JFunction1D_t::argument_type>                 JArray_t;

  JPDF4d_t i_pdf;

  NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
  {
    const string file_name = getFilename(fileDescriptor, SCATTERED_LIGHT_FROM_EMSHOWER);

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

    i_pdf.load(file_name.c_str());

    NOTICE("OK" << endl);
  }
  NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
  {
    JPDF4d_t pdf;

    const string file_name = getFilename(fileDescriptor, DIRECT_LIGHT_FROM_EMSHOWER);

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

    pdf.load(file_name.c_str());

    NOTICE("OK" << endl);
    NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);

    pdf.setExceptionHandler(new JPDF4d_t::JDefaultResult(JMATH::zero));

    NOTICE("Add PDF... " << flush);

    JTimer timer;

    timer.start();

    i_pdf.add(pdf);

    timer.stop();

    if (debug >= debug_t) {
      timer.print(cout, unit_t);
    }

    NOTICE("OK" << endl);
  }
  NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);

  i_pdf.setExceptionHandler(new JPDF4d_t::JDefaultResult(JMATH::zero));

  if ( TTS > 0.0 ){

    NOTICE("Smearing combined table with sigma [ns] = " << TTS << " ..." << flush);

    i_pdf.blur( TTS );

    NOTICE("OK" << endl);
  }

  JPDF5d_t o_pdf;
  
  const JFunction4DTransformer_t* i_transformer = dynamic_cast<JFunction4DTransformer_t*>(i_pdf.transformer->clone());
  const JFunction5DTransformer_t* o_transformer = (i_transformer != NULL ? new JFunction5DTransformer_t(*i_transformer) : NULL);

  if (o_transformer != NULL && debug >= debug_t) {
    o_transformer->print(cout);
  }

  if (Es.empty()) {
    Es.insert( 1.0 );
    Es.insert( 2.0 );
    Es.insert( 3.0 );
    Es.insert( 3.5 );
    Es.insert( 4.0 );
    Es.insert( 4.5 );
    Es.insert( 5.0 );
    Es.insert( 5.5 );
    Es.insert( 6.0 );
    Es.insert( 6.5 );
    Es.insert( 7.0 );
    Es.insert( 7.5 );
    Es.insert( 8.0 );
  };

  if (Ds.empty()) {
    Ds.insert(  0.10);
    Ds.insert(  0.50);
    Ds.insert(  1.00);
    Ds.insert(  2.00);
    Ds.insert(  5.00);
    Ds.insert(  8.00);
    Ds.insert( 10.00);
    Ds.insert( 20.00);
    Ds.insert( 30.00);
    Ds.insert( 40.00);
    Ds.insert( 50.00);
    Ds.insert( 60.00);
    Ds.insert( 70.00);
    Ds.insert( 80.00);
    Ds.insert( 90.00);
    Ds.insert(100.00);
    Ds.insert(120.00);
    Ds.insert(150.00);
    Ds.insert(170.00);
    Ds.insert(190.00);
    Ds.insert(210.00);
    Ds.insert(230.00);
    Ds.insert(250.00);
    Ds.insert(270.00);
    Ds.insert(290.00);
    Ds.insert(310.00);
    Ds.insert(400.00);
    Ds.insert(800.00);   // ARCA is big, we go to 800 m
  }

  set<double> C;        // cosine emission angle

  JQuadrature qeant(-1.0, +1.0, cosine_angle_count, geanx);

  for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {
    C.insert(i->getX());
  }
  
  C.insert(-1.00);
  C.insert(+1.00);
  
  set<double> X;        // time [ns]
  
  X.insert(  0.00);
  X.insert(  0.50);
  X.insert(  1.00);
  X.insert(  1.50);
  X.insert(  2.00);
  X.insert(  2.50);
  X.insert(  3.00);
  X.insert(  3.50);
  X.insert(  4.00);
  X.insert(  5.0);
  X.insert(  6.0);
  X.insert(  8.0);
  X.insert( 10.0);
  X.insert( 15.0);
  X.insert( 20.0);
  X.insert( 25.0);
  X.insert( 30.0);
  X.insert( 40.0);
  X.insert( 50.0);
  X.insert( 60.0);
  X.insert( 80.0);
  X.insert(100.0);
  X.insert(120.0);
  X.insert(150.0);
  X.insert(200.0);
  X.insert(250.0);
  X.insert(300.0);
  X.insert(350.0);
  X.insert(400.0);
  X.insert(450.0);
  X.insert(500.0);
  X.insert(600.0);
  X.insert(700.0);
  X.insert(800.0);
  X.insert(900.0);
  X.insert(1200.0);
  X.insert(1500.0);
    
  const double grid  = 7.0;                                       // [deg]  
  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0));  // azimuth angle unit step size
  
  NOTICE("building multi-dimensional function object" << endl);
      
  for (const double E : Es) {

    DEBUG("E: " << E << endl);
	
    vector< double > elong_lengths;

    for (size_t i = 0; i != elong_sample_count ; ++i ){

      const double fraction = (double) (i + 0.5) / (double) elong_sample_count;
      const double len      = geanz.getLength(pow(10.0, E), fraction, 1e-6);

      DEBUG("fraction: " << fraction << ", length [m]: "<< len << endl);

      elong_lengths.push_back(len);
    }

    if (elong_lengths.empty()) {
      elong_lengths.push_back(0.0);
    }

      
    for (const double D_m : Ds) {

      DEBUG("D_m: " << D_m << endl);
      
      for (const double cd : C) {

	NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);
	
	const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));
	
	for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {
	 
	  const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
	  
	  for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {
		
	    JFunction1D_t& f1 = o_pdf[E][D_m][cd][theta][phi];

	    f1.reserve(X.size());
	    
	    const JArray_t array(E, D_m, cd, theta, phi);

	    double t_old = (o_transformer != NULL ? o_transformer->getXn(array, *X.begin()) : *X.begin());
	    double y_old = 0.0;

	    for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
	  
	      const double t = (o_transformer != NULL ? o_transformer->getXn(array, *x) : *x);

	      double y = 0.0;
	      
	      for (const double Z : elong_lengths) {
	
		const double Z0 = cd * D_m;
		const double Z1 = Z0 - Z;
		
		const double D_m_prime = sqrt(D_m * D_m + Z1 * Z1 - Z0 * Z0);
		const double cd_prime  = Z1 / D_m_prime;
		const double t_prime   = t - (Z) / getSpeedOfLight() - (D_m_prime - D_m) * getIndexOfRefraction() / getSpeedOfLight();
		
		y += i_pdf(D_m_prime, cd_prime, theta, phi, t_prime);
	      }

	      y /= (double) elong_lengths.size();
	      
	      if (y != 0.0) {
		
		if (*x < 0.0) {
		  WARNING("dt < 0 " << *x << ' ' << D_m << ' ' << t << ' ' << y << endl);
		}
		
		if (y_old == 0.0) {
		  f1[t_old] = y_old;
		}
	      
		f1[t] = y;
		
	      } else {
		
		if (y_old != 0.0) {
		  f1[t] = y;
		}
	      }
	    
	      t_old = t;
	      y_old = y;
	    }
	  }
	}
      }
    }
  }
  NOTICE("OK" << endl);
  NOTICE("Memory " << FIXED(5,1) << getMemoryUsage() << "%" << endl);

  o_pdf.compile();

  if (o_transformer != NULL && option) {

    NOTICE("transform... " << flush);

    o_pdf.transform(*o_transformer);

    NOTICE("OK" << endl);
  }

  try {

    NOTICE("storing output to file " << outputFile << "... " << flush);

    o_pdf.store(outputFile.c_str());

    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
}