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

#include "TRandom3.h"
#include "TMath.h"

#include "JSirene/JSireneToolkit.hh"

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


/**
 * \file
 *
 * Example program to verify JSireneToolkit.hh.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  double   precision;
  double   probability;
  UInt_t   seed;
  int      debug;

  try { 

    JParser<> zap("Example program to verify JSireneToolkit.hh.");
    
    zap['e'] = make_field(precision)           = 1.0e-5;
    zap['P'] = make_field(probability)         = 1.0e-4;
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 3;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  gRandom->SetSeed(seed);

  {
    const size_t M = 100000;

    for (const double NPE : { 1.0e-8, 1.0e-5, 1.0e-3, 1.0e-2, 1.0e-1, 1.0, 1.0e+1, 1.0e+2 }) {

      const size_t numberOfEvents = max((size_t) (10.0 / NPE), M);

      DEBUG("NPE " << SCIENTIFIC(12,3) << NPE << ' ' << setw(10) << numberOfEvents << endl);

      vector<double> npe = { 0.001, 0.005, 0.010, 0.020, 0.050, 0.100, 0.200, 0.300 };

      double factor = 0.0;

      for (vector<double>::iterator i = npe.begin(); i != npe.end(); ++i) {
	factor += *i;
	*i     *= NPE;
      }
      
      ASSERT(factor < 1.0, "Test of relative number of photo-electrons " << FIXED(7,3) << factor);

      const size_t number_of_pmts = npe.size();

      vector<int> V(number_of_pmts, 0);

      for (size_t counter = 0; counter != numberOfEvents; ++counter) {

	if (counter%M == 0) {
	  STATUS(setw(10) << counter << '\r');
	}

	const size_t N = gRandom->Poisson(NPE);

	if (N != 0) {
	  
	  const vector<size_t>& ns = getNumberOfPhotoElectrons(NPE, N, npe, NPE < 25.0);
	  
	  size_t W = 0;

	  for (size_t i = 0; i != number_of_pmts; ++i) {
	    V[i] += ns[i];
	    W    += ns[i];
	  }

	  if (NPE < 25.0 && W > N) {
	    FATAL("Test of number of photo-electrons failed " << W << " > " << N << endl);
	  }
	}
      }
      STATUS(endl);

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

	const double x = (double) V[i]  / (double) numberOfEvents;

	DEBUG("npe " << SCIENTIFIC(12,3) << npe[i] << ' ' << SCIENTIFIC(12,3) << x << ' ' << SCIENTIFIC(12,3) << TMath::Poisson(x,npe[i]) << endl);

	ASSERT(TMath::Poisson(x,npe[i]) >= probability, "Test number of photo-electrons " << i);
      }
    }
  }
  {
    const double E  =  1.0e3;

    const double Tx[] = { +1.0e-3, -1.0e-3 };
    const double Ty[] = { -1.0e-3, +1.0e-3 };

    const double step = 10.0;

    JVertex vertex(0.0, 0.0, E);

    JTrack  muon(vertex);

    for (int i = 0; i != 2; ++i) {

      vertex.applyEloss(JVersor3Z(Tx[i],Ty[i]), 0.0);

      vertex.step(0.0, step);

      DEBUG("vertex "
	    << FIXED(7,3) << vertex.getX()  << ' '
	    << FIXED(7,3) << vertex.getY()  << ' '
	    << FIXED(7,3) << vertex.getZ()  << ' '
	    << FIXED(7,4) << vertex.getDX() << ' '
	    << FIXED(7,4) << vertex.getDY() << endl);

      muon.push_back(vertex);
    }

    {
      JVector2D pos = muon.getPosition(0.5*step);

      const double x = 0.0 + 0.5*step*Tx[0];
      const double y = 0.0 + 0.5*step*Ty[0];

      ASSERT(fabs(pos.getX() - x) <= precision);
      ASSERT(fabs(pos.getY() - y) <= precision);
    }

    {
      JVector2D pos = muon.getPosition(1.5*step);

      const double x = 0.0 + 1.0*step*Tx[0] + 0.5*step*(Tx[0] + Tx[1]);
      const double y = 0.0 + 1.0*step*Ty[0] + 0.5*step*(Ty[0] + Ty[1]);

      ASSERT(fabs(pos.getX() - x) <= precision);
      ASSERT(fabs(pos.getY() - y) <= precision);
    }
  }

  return 0;
}