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

#include "TROOT.h"
#include "TFile.h"

#include "jeep/parser.hh"
#include "jeep/timing.hh"
#include "antcc/Header.hh"
#include "antcc/Event.hh"
#include "trigger/trigger_io_helper.hh"
#include "trigger/algorithm.hh"
#include "trigger/cluster.hh"
#include "dataformat/AntaresTree.hh"

#include "geometry.hh"

#include "eureka.hh"


/**
 * Auxiliary class to sort hits in time (earliest hit first)
 */ 
class hit_sorter {
public:
  hit_sorter()
  {}

  bool operator()(const GEOMETRY::hit& x, const GEOMETRY::hit& y) const
  {
    return x.t() < y.t();
  }
};


/**
 * Auxiliary method to calculate chi2 for a giceb track and set of hits
 */
template<class T>
static double chi2(const GEOMETRY::track1Z& ta, T __begin, T __end)
{
  double chi2_ = 0;

  for (T i = __begin; i != __end; ++i) {
    
    double dx = i->x() - ta.x();
    double dy = i->y() - ta.y();
    double dz = i->z() - ta.z();
    double r  = sqrt(dx*dx + dy*dy);
    double dt = i->t() - ta.t() - (dz + r * UTIL::getTanThetaC()) * UTIL::C_INVERSE;
    
    chi2_ += (dt*dt) / (i->sigma()*i->sigma());
  }

  return chi2_;
}


/**
 * hit with associated address (and integer for fast permutations)
 */
class c_hit : public GEOMETRY::hit, public TRIGGER::integer, public TRIGGER::ars_address {
public:
  c_hit() :
    hit(),
    integer(),
    ars_address()
  {}

  c_hit(const GEOMETRY::hit& x, const TRIGGER::ars_address& id) :
    hit(x),
    integer(),
    ars_address(id)
  {}

  c_hit& operator=(const int& v_)
  {
    integer::operator=(v_);
    return *this;
  }
};


/**
 * Axiliary method to compare addresses
 */
bool address_comparator(const c_hit& x, const c_hit& y)
{
  return (const TRIGGER::ars_address&) x == (const TRIGGER::ars_address&) y;
}


/**
 * Auxiliary class to apply 1D match criterion after rotation
 */
class match1R : private UTIL::rotation, private CLUSTER::match1D {
public:
  match1R(const UTIL::direction& dir) :
    rotation(dir),
    match1D()
  {}

  template<class T> bool operator()(const T& _a, const T& _b)
  {
    T a(_a);
    T b(_b);

    a.rotate(*this);
    b.rotate(*this);

    return CLUSTER::match1D::operator()(a,b);
  }
};


/**
 * program to make linear prefits in various directions
 * input   ROOT  Event
 *               PhysicsEvent
 * output  ASCII Event + prefit results
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JEEP;
  using namespace EUREKA;
  using TRIGGER::ars_address;
  using TRIGGER_IO::SPE_reader;
  using TRIGGER_IO::TriggerInterfaceHelper;
  using GEOMETRY::track1Z;
  using GEOMETRY::hit;

  string         inputFile;
  string         outputFile;
  string         interfaceInput;
  unsigned int   numberOfEvents;
  float          gridAngle;
  unsigned int   numberOfRemovableHits;
  bool           addL0;
  bool           debug;

  try {

    Parser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)            = "scanf.dat";
    zap['a'] = make_field(interfaceInput);
    zap['n'] = make_field(numberOfEvents)        = numeric_limits<int>::max();
    zap['g'] = make_field(gridAngle)             = 10.0;    // [deg]
    zap['r'] = make_field(numberOfRemovableHits) = 3; 
    zap['A'] = make_field(addL0);
    zap['d'] = make_field(debug);

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


  TriggerInterfaceHelper interface(interfaceInput.c_str());

  const float Max_Transverse_Distance = CLUSTER::getMaxTransverseDistance();

  Omega omega(gridAngle);

  
  ofstream out(outputFile.c_str());


  {
    TFile* in = TFile::Open(inputFile.c_str());
    
    Header* header = (Header*) in->Get(Header::Class_Name());
    
    interface.offset(position(-header->coord_origin.x,
			      -header->coord_origin.y,
			      -header->coord_origin.z));
    
    in->Close();
  }


  MonteCarloEvent_Reader.Add(inputFile.c_str());
  PhysicsEvent_Reader.Add(inputFile.c_str());


  const double        SIGMA                   =  5; // [ns]
  const unsigned int  MINIMUM_NUMBER_OF_HITS  =  5;


  size_t physicsEventCounter   = 0;
  size_t numberOfPhysicsEvents = (size_t) PhysicsEvent_Reader.GetEntries();

  PhysicsEvent* result = NULL;


  numberOfEvents = min(numberOfEvents, (unsigned int) MonteCarloEvent_Reader.GetEntries());

  for (unsigned int eventCounter = 0; eventCounter < numberOfEvents; ++eventCounter) {

    MonteCarloEvent_Reader.GetEvent(eventCounter);

    Event* event = MonteCarloEvent_Reader.GetAddress();

    cerr << "event: " << setw(8) << event->eventNumber() << '\r';


    // Find physics event (check!)

    for (; (result == NULL || (int) result->TriggerCounter() < event->eventNumber()) &&
           physicsEventCounter < numberOfPhysicsEvents;
         ++physicsEventCounter) {

      PhysicsEvent_Reader.GetEvent(physicsEventCounter);

      result = PhysicsEvent_Reader.GetAddress();
    }


    if (result != NULL && (int) result->TriggerCounter() == event->eventNumber()) {


      // ASCII output
      
      Event::OStream os(out,*event);
      
      
      double timeOfRTS = CLOCK::getTimeOfRTS(event->eventTime());

      vector<c_hit> data_L1;
      vector<c_hit> data_L0;
    

      for (PhysicsEvent::const_iterator<TriggeredSPE_Hit> i = result->begin<TriggeredSPE_Hit>(); i != result->end<TriggeredSPE_Hit>(); ++i) {

	ars_address id(i->lcm_id, i->ars_id);

	const SPE_reader& spe_reader = interface.getSPEreader(id);

	double t1 = timeOfRTS  +  spe_reader.timeSinceRTS(i->timestamp,i->tvc);
	
	hit o(interface.getPosition(id), t1, SIGMA);

	data_L1.push_back(c_hit(o,id));
      }


      for (PhysicsEvent::const_iterator<SPE_Hit> i = result->begin<SPE_Hit>(); i != result->end<SPE_Hit>(); ++i) {

	ars_address id(i->lcm_id, i->ars_id);

	const SPE_reader& spe_reader = interface.getSPEreader(id);

	double t1 = timeOfRTS  +  spe_reader.timeSinceRTS(i->timestamp,i->tvc);
	
	hit o(interface.getPosition(id), t1, SIGMA);

	data_L0.push_back(c_hit(o,id));
      }

    
      s_timing timer;
    
      timer.reset(); 
      timer.start(); 
    
    
      for (t_omega::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
      
      
	CLUSTER::setMaxTransverseDistance(Max_Transverse_Distance * (1.0 - 0.25*dir->cost()*dir->cost()));
      

	vector<c_hit>::const_iterator __begin = data_L1.begin();
	vector<c_hit>::const_iterator __end   = CLUSTER::cluster1D(data_L1.begin(), data_L1.end(), *dir);


	vector<c_hit> buffer;
	  

	if (addL0) {

	  match1R match(*dir);
	  
	  for (vector<c_hit>::const_iterator i = data_L0.begin(); i != data_L0.end(); ++i) {
	    
	    bool ok = true;
	    
	    for (vector<c_hit>::const_iterator j = __begin; j != __end && ok; ++j)
	      ok = !address_comparator(*i,*j) && match(*i,*j);
	    
	    if (ok)
	      buffer.push_back(*i);
	  }


	  buffer.erase(CLUSTER::cluster1D(buffer.begin(), buffer.end(), *dir), buffer.end());
	

	  for (vector<c_hit>::const_iterator i = __begin; i != __end; ++i)
	    buffer.push_back(*i);

	} else {

	  copy(__begin, __end, back_inserter(buffer));
	}
 

	// earliest hit in each storey

	if (true) {

	  map<int, set<c_hit,hit_sorter> > wmap;
	
	  for (vector<c_hit>::iterator i = buffer.begin(); i != buffer.end(); ++i)
	    wmap[i->lcm_id()].insert(*i);

	  buffer.clear();
	
	  for (map<int, set<c_hit,hit_sorter> >::const_iterator i = wmap.begin(); i != wmap.end(); ++i)
	    buffer.push_back(*(i->second.begin()));
	}

      
	if (buffer.size() >= MINIMUM_NUMBER_OF_HITS) {
	
	  // rotate
	
	  rotation R(*dir);

	  for (vector<c_hit>::iterator i = buffer.begin(); i != buffer.end(); ++i)
	    i->rotate(R);
	
	  try {
	  
	    track1Z tb;
	  
	    using namespace TRIGGER;
	  
	    generate(buffer.begin(), buffer.end(), integer(0));
	  
	    const unsigned int n = min(buffer.size() - MINIMUM_NUMBER_OF_HITS, numberOfRemovableHits);
	  
	    vector<c_hit>::iterator q;
	  
	    advance(q = buffer.begin(), buffer.size() - n);
	  
	    double chi2_min = numeric_limits<double>::max();
	  
	    do {
	    
	      try {
	      
		track1Z tc    = eureka(buffer.begin(), q);
		double  chi2_ = chi2(tc, buffer.begin() ,q);
	      
		if (chi2_ < chi2_min) {
		  tb       = tc;
		  chi2_min = chi2_;
		}
	      }
	      catch(exception& error) { }
	    
	    } while (next_permutation(buffer.begin(), q, buffer.end()));

	    if (chi2_min <  numeric_limits<double>::max())
	      os << "eureka: " << tb << ' ' << *dir << ' ' << chi2_min << ' ' << buffer.size() - n << endl;
	  }
	  catch(exception& error) {
	    cerr << endl << error.what() << endl;
	  }
	}
      }

      timer.stop(); 

      os << "timing: " << (double) timer.usec_ucpu << endl;
    }
  }
  cerr << endl;

  out.close();
}