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

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

#include "jeep/parser.hh"
#include "antcc/Header.hh"
#include "antcc/Event.hh"
#include "trigger/util.hh"

#include "geometry.hh"
#include "tools.hh"


/**
 * data structure for eureka results (output of scanf)
 */
class eureka : public GEOMETRY::track1Z, public UTIL::direction {
public:
  eureka()
  {}
  
  eureka(std::istream& in)
  {
    read(in);
  }
    
  std::istream& read(std::istream& in)
  {
    GEOMETRY::track1Z::read(in);
    UTIL::direction::read(in);
    in >> chi2_ >> N_;
    return in;
  }

  const double& chi2() const { return chi2_; }
  const int&    N()    const { return N_; }

protected:
  double chi2_;
  int    N_;
};


/**
 * tentative sorter of solutions; to be optimised?
 */
static inline bool eureka_sorter(const eureka& x, const eureka& y)
{
  return (x.N() - x.chi2()/x.N())  >  (y.N() - y.chi2()/y.N());
  //return x.chi2()/x.N()  <  y.chi2()/y.N();
}


/**
 * Extension of Event class to include eureka results
 */
class ZEvent : public Event {
public:
  ZEvent() :
    Event()
  {}

  ZEvent(std::istream& in) :
    Event()
  {
    Streamer(in);
  }

  virtual bool filter(std::istream& in, const std::string& key, const bool& isReadingEvent) 
  { 
    if (isReadingEvent) {

      if      (key == "eureka:")
	result.push_back(eureka(in));
      else if (key == "timing:")
	in >> usec;
      else
	return false; 
    }
     
    return true;
  }

  void clear()
  {
    Event::clear();
    result.clear();
    usec = 0;
  }

  std::vector<eureka> result;
  double              usec;
};


/**
 * Analysis program for scanf results
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace UTIL;
  using namespace TOOLS;
  using namespace GEOMETRY;

  string         inputFile;
  string         outputFile;
  unsigned int   numberOfEvents;
  unsigned int   numberOfSolutions;
  bool           debug;

  try {

    Parser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)        = "doemaar.root";
    zap['n'] = make_field(numberOfEvents)    = numeric_limits<int>::max();
    zap['N'] = make_field(numberOfSolutions) = 1;
    zap['d'] = make_field(debug);

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


  TFile out(outputFile.c_str(), "recreate");

  TH1F hn("n", "index        ",  21, -0.5,  20.5);
  TH1F ha("a", "delta a [deg]", 500,  0.1, 180.0);
  TH1F hp("p", "probability  ", 100,  0.0,   1.0);
  TH1F hu("u", "cpu time [us]", 100,  0.0,   1.0e6);

  ifstream in(inputFile.c_str());

  unsigned int n = 0;

  for (ZEvent event; event.Streamer(in); event.clear()) {

    cerr << "event: " << event.eventNumber() << '\r';

    if (++n > numberOfEvents)
      break;

    // Monte Carlo muon

    vector<McTrack>::const_iterator p;
    
    for (p = event.TrackList().begin(); p != event.TrackList().end(); ++p) {
      if (p->type() == 5)
        break;
    }

    if (p != event.TrackList().end() && event.result.size() != 0) {


      direction dir(p->direction().x, p->direction().y, p->direction().z);


      // sort; best solution first

      sort(event.result.begin(), event.result.end(), eureka_sorter);


      // find closest solution within the first few

      double dmin = -1;
      int    j    =  0;

      for (unsigned int i = 0; i != min(numberOfSolutions,event.result.size()); ++i) {

	const eureka& ta = event.result[i];

	double d = dot(dir,ta);
	  
	if (abs(1.0 - d) < abs(1.0 - dmin)) {
	  j    = i;
	  dmin = d;
	}
      }

      // histogram

      const eureka& ta = event.result[j];
	  
      double d = dot(dir,ta);
      double p = NR::prob(ta.chi2(), ta.N() - 3);

      hn.Fill((double) j);
      ha.Fill(acos(d) * 180 / PI);
      hp.Fill(p);
      hu.Fill(event.usec);
    }
  }
  cerr << endl;

  in.close();

  out.Write();
  out.Close();
}