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

#include "TRandom3.h"

#include "JTools/JQuantile.hh"
#include "JTrigger/JMatch3B.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JGeometry3D/JSphere3D.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JVertex3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JGeometry3DTestkit.hh"

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


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

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

  try { 

    JParser<> zap("Example program to test hit clusterize with JMatch3B.");
    
    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['P'] = make_field(P)                   =   0.1;
    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 JMatch3B<JHit_t>    match(roadWidth_m, Tmax_ns);

  {
    // signal

    const JSphere3D sphere(JVector3D(0,0,0), 0.5*roadWidth_m);
    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), getRandom<JDirection3D>(), 0.0);

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

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

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

      ASSERT(clusterize(data.begin(), data.end(), match) == data.end());
    }
  }

  {
    // background

    const JSphere3D sphere(JVector3D(0,0,0), 5.0*roadWidth_m);
    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(sphere);
	const double    t1  = gRandom->Uniform(-T1max, +T1max);

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

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

    DEBUG(Q);

    ASSERT(Q.getMean() < P);
  }

  return 0;
}