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

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

#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"

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


/**
 * \file
 *
 * Example program to histogram time-over-threshold distributions.
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JPMTAnalogueSignalProcessor pmt;
  string                      outputFile;
  double                      mu;
  int                         debug;

  try {

    JProperties properties = pmt.getProperties();

    JParser<> zap("Example program to histogram time-over-threshold distributions.");

    zap['o'] = make_field(outputFile)          = "pmt.root";
    zap['P'] = make_field(properties)          = JPARSER::initialised();
    zap['E'] = make_field(mu);
    zap['d'] = make_field(debug)               = 2;

    if (zap.read(argc, argv) != 0)
      return 1;
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  pmt.configure();

  NOTICE("Time-over-threshold [ns] " << pmt.getDecayTime(1.0, pmt.threshold) + pmt.getRiseTime(1.0, pmt.threshold) << endl);

  const double x1 = pmt.getStartOfLinearisation();
  const double z  = pmt.riseTime_ns / sqrt(2.0 * log(x1/pmt.threshold));

  NOTICE("Transition point [npe]   = " << FIXED(6,3) << x1                            << endl);
  NOTICE("Transition point [ns]    = " << FIXED(6,3) << pmt.getTimeOverThreshold(x1)  << endl);
  NOTICE("d(npe)/d(T) (exact)      = " << FIXED(6,4) << 1.0 / pmt.slope               << endl); 
  NOTICE("d(npe)/d(T) (calculated) = " << FIXED(6,4) << x1 / (z + pmt.getDecayTime()) << endl);


  TH1D h0("h0", NULL, 1200, -30.0, +60.0);

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

    const double t1_ns = h0.GetBinCenter(i);
    
    h0.SetBinContent(i, pmt.getAmplitude(t1_ns));
  }

  const Double_t X = pmt.getT1();
  const Double_t Y = pmt.getY1();

  DEBUG("Transition point (" << FIXED(5,2) << X << "," << FIXED(5,2) << Y << ")" << endl);

  TGraph g0(1, &X, &Y);

  g0.SetName("g0");

  Double_t x[4];
  Double_t y[4];

  x[0] = -pmt.getRiseTime(1.0, pmt.threshold);
  y[0] =  0.0;
  x[1] = -pmt.getRiseTime(1.0, pmt.threshold);
  y[1] =  pmt.threshold;
  x[2] =  pmt.getDecayTime(1.0, pmt.threshold);
  y[2] =  pmt.threshold;
  x[3] =  pmt.getDecayTime(1.0, pmt.threshold);
  y[3] =  0.0;
  
  TGraph g1(4, x, y);

  g1.SetName("g1");

  JManager<int, TH1D> zmap(new TH1D("h1[%]", NULL, 4000, 0.0, +70.0));

  for (int NPE = 1; NPE != 7; ++NPE) {

    for (Int_t i = 1; i <= zmap->GetNbinsX(); ++i) {

      const double tot_ns = zmap->GetBinCenter(i);
      const double npe    = pmt.getNPE(tot_ns);

      const double W      = TMath::PoissonI((Double_t) NPE, mu) * pmt.getChargeProbability(npe, NPE) * pmt.getDerivative(npe);

      zmap[ 0 ]->AddBinContent(i, W);
      zmap[NPE]->SetBinContent(i, W);
    }
  }


  TH1D h2("h2", NULL, 5000, 0.0, 1000.0);

  for (Int_t i = 1; i <= h2.GetNbinsX(); ++i) {

    const double npe = h2.GetBinCenter(i);

    h2.SetBinContent(i, pmt.getTimeOverThreshold(npe));
  }


  TH1D h3("h3", NULL, 5000, 0.0, +5.0);

  const int NPE = 1;

  for (Int_t i = 1; i <= h3.GetNbinsX(); ++i) {

    const double x = h3.GetBinCenter(i);
    const double y = pmt.getChargeProbability(x, NPE);

    h3.SetBinContent(i, y);
  }


  TH1D h4("h4", NULL, 5000, 0.0, +5.0);

  for (Int_t i = 1; i <= h4.GetNbinsX(); ++i) {

    const double x = h4.GetBinCenter(i);
    const double y = pmt.getDerivative(x);

    h4.SetBinContent(i, y);
  }


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

  out << h0 << g0 << g1 << zmap << h2 << h3 << h4;

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