#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>

#include "TRandom3.h"

#include "JTools/JQuantile.hh"
#include "JTools/JRange.hh"
#include "JTrigger/JMatch1D.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JVertex3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JGeometry3DTestkit.hh"

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


/**
 * \file
 * Example program to test hit JTRIGGER::clusterize with JTRIGGER::JMatch1D.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef  JRange<double>  JRange_t;

  int      numberOfEvents;
  int      numberOfPoints;
  double   Tmax_ns;
  double   roadWidth_m;
  JRange_t Z_m;
  double   P;
  UInt_t   seed;
  int      debug;

  try { 

    JParser<> zap("Example program to test hit clusterize with JMatch1D.");
    
    zap['n'] = make_field(numberOfEvents)      = 100; 
    zap['N'] = make_field(numberOfPoints)      = 100; 
    zap['T'] = make_field(Tmax_ns)             =  15.0;   // [ns]
    zap['R'] = make_field(roadWidth_m)         = 100.0;   // [m]
    zap['z'] = make_field(Z_m)                 = JRange_t(-500.0, +500.0);  // [m]
    zap['P'] = make_field(P)                   =   0.05;
    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);


  typedef JVertex3D         JHit_t;
  typedef vector<JHit_t>    JData_t;

  const JMatch1D<JHit_t>    match(roadWidth_m, Tmax_ns);

  {
    // signal

    JTimer timer("Signal");

    const JCylinder3D cylinder(JCircle2D(JVector2D(0,0),0.5*roadWidth_m), Z_m.getLowerLimit(), Z_m.getUpperLimit());
    const double      T1max = 0.45*Tmax_ns;

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

      STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);

      const JTrack3D trk(JVector3D(0,0,0), JVersor3D(0,0,1), 0.0);

      JData_t data;
    
      for (int i = 0; i != numberOfPoints; ++i) {

	const JVector3D pos = getRandomPosition(cylinder);
	const double    t1  = trk.getT(pos) + gRandom->Uniform(-T1max, +T1max);

	data.push_back(JHit_t(pos, t1));
      }

      timer.start();

      JData_t::iterator __end = clusterize(data.begin(), data.end(), match);

      timer.stop();

      //ASSERT(__end == data.end(), "Test of signal.");
    }

    timer.print(cout, true, micro_t);
  }

  {
    // background

    JTimer timer("Background");

    const JCylinder3D cylinder(JCircle2D(JVector2D(0,0),5.0*roadWidth_m), Z_m.getLowerLimit(), Z_m.getUpperLimit());
    const double      T1max = 5.0*roadWidth_m / C + Tmax_ns;

    JQuantile Q("background");

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

      STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);

      JData_t data;

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

	const JVector3D pos = getRandomPosition(cylinder);
	const double    t1  = gRandom->Uniform(-T1max, +T1max);

	data.push_back(JHit_t(pos, t1));
      }

      timer.start();

      JData_t::iterator __end = clusterize(data.begin(), data.end(), match);

      timer.stop();

      Q.put((double) distance(data.begin(), __end) / (double) data.size());
    }
    STATUS(endl);

    timer.print(cout, true, micro_t);

    DEBUG(Q);

    ASSERT(Q.getMean() < P, "Test of background.");
  }

  return 0;
}