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

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

#include "JPhysics/JPDF.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JConstants.hh"

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

#include <cmath>

#include "JTools/JFunction1D_t.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/JPDFToolkit.hh"

namespace {
 
  using JTOOLS::JPolint1Function1D_t;
 
  /**
   * Auxiliary class for correction of number of photons from EM-shower.
   */
  class JEMShowerCorrection
  {
  public:
    /**
     * Default constructor.
     */
    JEMShowerCorrection() 
    {
      using namespace JPP;

      const double P_Atm = 240.0;                      // ambient pressure [Atm]
      const double wmin  = getMinimalWavelength();     // [nm]
      const double wmax  = getMaximalWavelength();     // [nm]

      const JDispersion dispersion(P_Atm);

      const double xmin = 1.0 / wmax;
      const double xmax = 1.0 / wmin;
      const int    nx   = 10000;

      double value = 0.0;

      for (double x = xmin, dx = (xmax - xmin) / (double) nx; x <= xmax; x += dx) {

        const double w  = 1.0 / x;
        const double dw = dx * w*w;

        const double n  = dispersion.getIndexOfRefractionPhase(w);
        
        value += cherenkov(w,n) * dw;
      }

      value *= geanc();                                // number of photons per GeV


      f1[MASS_ELECTRON]  =  0.00;
      f1[ 1.0e-3]  =    90.96 /  1.0e-3;               // number of photons per GeV
      f1[ 2.0e-3]  =   277.36 /  2.0e-3;               // from Geant4 simulation by Daniel Lopez
      f1[ 3.0e-3]  =   485.82 /  3.0e-3;
      f1[ 4.0e-3]  =   692.83 /  4.0e-3;
      f1[ 5.0e-3]  =   890.01 /  5.0e-3;
      f1[ 6.0e-3]  =  1098.53 /  6.0e-3;
      f1[ 7.0e-3]  =  1285.47 /  7.0e-3;
      f1[ 8.0e-3]  =  1502.86 /  8.0e-3;
      f1[ 9.0e-3]  =  1687.15 /  9.0e-3;
      f1[10.0e-3]  =  1891.00 / 10.0e-3;

      f1.div(value);
      f1.compile();

      f2[-2.0]  =  log(1.891e3 / 1.0e-2);
      f2[-1.0]  =  log(1.905e4 / 1.0e-1);
      f2[ 0.0]  =  log(1.889e5 / 1.0e+0);
      f2[+1.0]  =  log(1.875e6 / 1.0e+1);
      f2[+2.0]  =  log(1.881e7 / 1.0e+2);

      f2.sub(log(value));
      f2.compile();
    }


    /**
     * Get correction of number of photons from EM-shower as a function of energy.
     *
     * \param  E          energy [GeV]
     * \return            correction
     */
    double operator()(const double E) const
    {
      if        (E <= f1.getXmin()) {

	return 0.0;

      } else if (E <= f1.getXmax()) {

	return f1(E);

      } else {

	const double x = log10(E);

	if        (x <= f2.getXmin()) {

	  return exp(f2.begin()->getY());

	} else if (x <= f2.getXmax()) {

	  return exp(f2(x));

	} else {

	  return exp(f2.rbegin()->getY());
	}
      }
    }


  private:
    JPolint1Function1D_t f1;
    JPolint1Function1D_t f2;
  };
  

  /**
   * Function object for EM-shower correction
   */
  static const JEMShowerCorrection  getEMShowerCorrection;
}


/**
 * \file
 *
 * Auxiliary program to draw npe as a function of EM-energy.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;

  string      outputFile;
  int         numberOfPoints;
  double      epsilon;
  int         debug;

  try { 

    JParser<> zap("Auxiliary program to draw npe as a function of EM-energy.");

    zap['o'] = make_field(outputFile)      = "geanc.root";
    zap['n'] = make_field(numberOfPoints)  = 25;
    zap['e'] = make_field(epsilon)         = 1.0e-10;
    zap['d'] = make_field(debug)           = 2;

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

  const JPDF_C
    pdf(NAMESPACE::getPhotocathodeArea(),
	NAMESPACE::getQE,
	NAMESPACE::getAngularAcceptance,
	NAMESPACE::getAbsorptionLength,
	NAMESPACE::getScatteringLength,
        NAMESPACE::getScatteringProbability,
	NAMESPACE::getAmbientPressure(),
	getMinimalWavelength(),
	getMaximalWavelength(),
	numberOfPoints,
	epsilon);


  double xmin = -4.0;
  double xmax =  0.0;

  const double precision = 1.0e-10;

  while (fabs(xmax - xmin) > precision) {

    const double x = 0.5 * (xmin + xmax);
    const double E = pow(10.0, x);

    const double y = getEMShowerCorrection(E);

    if (y < 0.5) 
      xmin = x;
    else
      xmax = x;
  }

  const double EMIN =  pow(10.0, 0.5*(xmax + xmin));

  NOTICE("Threshold kinetic energy [GeV] " << sqrt((EMIN + MASS_ELECTRON) * (EMIN - MASS_ELECTRON)) << endl);

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

  TH1D h0("h0", NULL, 10000, -4.0, +4.0);
  TH1D h1("h1", NULL, 10000, -4.0, +4.0);
 
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {

    const Double_t x = h0.GetBinCenter(i);
    const Double_t E = pow(10.0, x);

    h0.SetBinContent(i, E * geanc() * pdf.getNumberOfPhotons());
    h1.SetBinContent(i, getEMShowerCorrection(E));
  }

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