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

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

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Auxiliary program to determine L0 and L1 hit probability as a function of 
 * the minimal distance of approach of a muon to a PMT.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string        fileDescriptor;
  string        outputFile;
  string        option;
  double        E_GeV;
  JAngle3D      dir;
  double        TMax_ns;
  int           debug;

  try { 

    JParser<> zap;

    zap['f'] = make_field(fileDescriptor);
    zap['o'] = make_field(outputFile)       = "multipmt.root";
    zap['O'] = make_field(option)           = "KM3NeT", "Antares";
    zap['E'] = make_field(E_GeV);
    zap['D'] = make_field(dir);
    zap['T'] = make_field(TMax_ns)          = 20.0;    // L1 coincidence time window [ns]
    zap['d'] = make_field(debug)            =  0;

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

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


  typedef JSplineFunction1S_t                                     JFunction1D_t;
  typedef
    JMapList<JPolint1FunctionalMap,
    JMapList<JPolint1FunctionalGridMap,
    JMapList<JPolint1FunctionalGridMap> > >                       JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  /**
   * PDF types.
   */
  const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
			       SCATTERED_LIGHT_FROM_MUON,
			       DIRECT_LIGHT_FROM_EMSHOWERS,
			       SCATTERED_LIGHT_FROM_EMSHOWERS };

  const double TTS            = 2.0;        // [ns]
  const int    numberOfPoints = 25;
  const double epsilon        = 1.0e-10;
 
  const int  NUMBER_OF_PDFS = sizeof(pdf_t)/sizeof(pdf_t[0]);

  JPDF_t pdf[NUMBER_OF_PDFS];               // PDF

  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));

  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
    
    try {
      
      const string file_name = getFilename(fileDescriptor, pdf_t[i]);
      
      NOTICE("loading PDF from file " << file_name << "... " << flush);
      
      pdf[i].load(file_name.c_str());
      
      NOTICE("OK" << endl);
      
      pdf[i].setExceptionHandler(supervisor);      
      pdf[i].blur(TTS, numberOfPoints, epsilon);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
  }


  // set-up

  vector<JAngle3D> PMT;

  if        (option == "KM3NeT") {

    PMT.push_back(JAngle3D(3.142, 0.000));   //  1
    PMT.push_back(JAngle3D(2.582, 0.524));
    PMT.push_back(JAngle3D(2.582, 1.571));
    PMT.push_back(JAngle3D(2.582, 2.618));
    PMT.push_back(JAngle3D(2.582, 3.665));
    PMT.push_back(JAngle3D(2.582, 4.712));
    PMT.push_back(JAngle3D(2.582, 5.760));
    PMT.push_back(JAngle3D(2.162, 0.000));
    PMT.push_back(JAngle3D(2.162, 1.047));
    PMT.push_back(JAngle3D(2.162, 2.094));   // 10
    PMT.push_back(JAngle3D(2.162, 3.142));
    PMT.push_back(JAngle3D(2.162, 4.189));
    PMT.push_back(JAngle3D(2.162, 5.236));
    PMT.push_back(JAngle3D(1.872, 0.524));
    PMT.push_back(JAngle3D(1.872, 1.571));
    PMT.push_back(JAngle3D(1.872, 2.618));
    PMT.push_back(JAngle3D(1.872, 3.665));
    PMT.push_back(JAngle3D(1.872, 4.712));
    PMT.push_back(JAngle3D(1.872, 5.760));
    PMT.push_back(JAngle3D(1.270, 0.000));   // 20
    PMT.push_back(JAngle3D(1.270, 1.047));
    PMT.push_back(JAngle3D(1.270, 2.094));
    PMT.push_back(JAngle3D(1.270, 3.142));
    PMT.push_back(JAngle3D(1.270, 4.189));
    PMT.push_back(JAngle3D(1.270, 5.236));
    PMT.push_back(JAngle3D(0.980, 0.524));
    PMT.push_back(JAngle3D(0.980, 1.571));
    PMT.push_back(JAngle3D(0.980, 2.618));
    PMT.push_back(JAngle3D(0.980, 3.665));
    PMT.push_back(JAngle3D(0.980, 4.712));   // 30
    PMT.push_back(JAngle3D(0.980, 5.760));

  } else if (option == "Antares") {

    PMT.push_back(JAngle3D(2.36, +1.05));
    PMT.push_back(JAngle3D(2.36,  3.14));
    PMT.push_back(JAngle3D(2.36, -1.05));

  } else {

    FATAL("Fatal error at detector option " << option << endl);
  };

    
  // rotate PMTs according specified orientation

  const JRotation3D R(dir);

  for (vector<JAngle3D>::iterator i = PMT.begin(); i != PMT.end(); ++i) {
    *i = JDirection3D(*i).rotate(R);
  }


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


  TH1D h0("L0", NULL, 300, 1.0, 151.0);
  TH1D h1("L1", NULL, 300, 1.0, 151.0);


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

    const double R = h0.GetBinCenter(i);

    double V = 0.0;

    for (vector<JAngle3D>::const_iterator pmt = PMT.begin(); pmt != PMT.end(); ++pmt) {

      V += (pdf[0](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V      +
	    pdf[1](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V      +
	    pdf[2](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V * E_GeV  +
	    pdf[3](R, pmt->getTheta(), fabs(pmt->getPhi()), 0.0).V * E_GeV);
    }

    h0.SetBinContent(i, 1.0 - TMath::PoissonI(0,V));
  }


  const int NUMBER_OF_PMTS = PMT.size();

  double Vi[NUMBER_OF_PMTS];     // integrals

  const double Tmin =  -15.0;    // [ns]
  const double Tmax = +250.0;    // [ns]
  const double dt   =    1.0;    // [ns]
    

  for (int i = 1; i <= h1.GetNbinsX(); ++i) {

    const double R = h1.GetBinCenter(i);

    double Y = 0.0;
      
    for (double x = Tmin; x <= Tmax; x += dt) {

      double V = 0.0;
      
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {

	Vi[pmt]  =  ((pdf[0](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v          +
		      pdf[1](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v          +
		      pdf[2](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v * E_GeV  +
		      pdf[3](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x + TMax_ns).v * E_GeV)
		     
		     -
		     
		     (pdf[0](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v          +
		      pdf[1](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v          +
		      pdf[2](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v * E_GeV  +
		      pdf[3](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).v * E_GeV));

	if (Vi[pmt] < 0.0) {
	  Vi[pmt] = 0.0;
	}
	
	V += Vi[pmt];
      }
      
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
	
	const double y = (pdf[0](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f          +
			  pdf[1](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f          +
			  pdf[2](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f * E_GeV  +
			  pdf[3](R, PMT[pmt].getTheta(), fabs(PMT[pmt].getPhi()), x).f * E_GeV);
	
	if (y > 0.0) {
	  Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
	}
      }  
    }
    
    h1.SetBinContent(i, 1.0 - TMath::PoissonI(0,Y));
  }

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