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

#include "TRandom3.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;
  size_t   numberOfEvents;
  double   sigma;
  UInt_t   seed;
  int      debug;

  try { 

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

  cout.tie(&cerr);

  gRandom->SetSeed(seed);

  ASSERT(numberOfEvents > 0);

  {
    const size_t number_of_pmts = 20;
    const double factor         = 0.5; // [0,1]

    for (const double NPE : { 0.2, 1.0, 10.0, 100.0 }) {

      const double npe = factor * NPE / number_of_pmts;

      double W[] = { 0.0, 0.0 };

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

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

	W[0] += N;

	if (N != 0) {
	  for (size_t i = 0; i != number_of_pmts; ++i) {
	    W[1] += getNumberOfPhotoElectrons(NPE, N, npe);
	  }
	}
      }

      W[0] /= numberOfEvents;
      W[1] /= numberOfEvents * number_of_pmts;

      DEBUG("NPE " << FIXED(7,3) << NPE << ' ' << FIXED(7,3) << W[0] << endl);
      DEBUG("npe " << FIXED(7,3) << npe << ' ' << FIXED(7,3) << W[1] << endl);

      ASSERT(fabs(NPE - W[0]) <= NPE * sigma / sqrt(numberOfEvents));
      ASSERT(fabs(npe - W[1]) <= npe * sigma / sqrt(numberOfEvents));
    }
  }
  {
    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;
}