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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"

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

#include "trigger/montecarlo_helper.hh"

#include "dataformat/AntaresTree.hh"


UTIL::position  getPosition (const Vec3D& v) { return UTIL::position (v.x, v.y, v.z); }
UTIL::direction getDirection(const Vec3D& v) { return UTIL::direction(v.x, v.y, v.z); }


/**
 * Auxiliary program to verify time residuals and PMT angular acceptance.
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace TRIGGER_IO;

  vector<string> inputFile;
  string         outputFile;
  string         interfaceInput;
  Long64_t       numberOfEvents;
  bool           debug;

  try { 

    Parser<> zap;
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "dot.root";
    zap['a'] = make_field(interfaceInput);
    zap['n'] = make_field(numberOfEvents)      = 1;
    zap['d'] = make_field(debug);
    
    if (zap.read(argc, argv) != 0)
      return 1;
  }
  catch(const exception &error) {
    cerr << error.what() << endl;
    return 2;
  }


  MonteCarloHelper interface(interfaceInput.c_str());

  if (debug)
    interface.print(cout);

  // ROOT

  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i)
    MonteCarloEvent_Reader.Add(i->c_str());

  Header header(MonteCarloEvent_Reader);

  if (header.getNumberOfHeaders() != 0) {

    Event::set_runNumber(header.start_run.run_id);

    interface.offset(position(-header.coord_origin.x,
                              -header.coord_origin.y,
                              -header.coord_origin.z));
  }




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


  map<int, TH1*> xmap;

  xmap[+5] = new TH1D("dt[+5]", NULL, 200, -50.0, +50.0);
  xmap[-5] = new TH1D("dt[-5]", NULL, 200, -50.0, +50.0);
  xmap[+3] = new TH1D("dt[+3]", NULL, 200, -50.0, +50.0);
  xmap[-3] = new TH1D("dt[-3]", NULL, 200, -50.0, +50.0);

  map<int, TH1*> ymap;

  ymap[+5] = new TProfile("ct[+5]", NULL, 100,  -1.0,  +1.0, -1e2, +1e2);
  ymap[-5] = new TProfile("ct[-5]", NULL, 100,  -1.0,  +1.0, -1e2, +1e2);
  ymap[+3] = new TProfile("ct[+3]", NULL, 100,  -1.0,  +1.0, -1e2, +1e2);
  ymap[-3] = new TProfile("ct[-3]", NULL, 100,  -1.0,  +1.0, -1e2, +1e2);

  TH1D h0("ct", NULL, 100,  -1.0,  +1.0);

  const double RMAX = 100.0;  // [m]
  const double RMIN =   5.0;  // [m]


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

  for (Long64_t event_count = 0; event_count < numberOfEvents; ++event_count) {

    MonteCarloEvent_Reader.GetEvent(event_count);

    Event* pevt = MonteCarloEvent_Reader.GetAddress();
    
    cerr << "event: " << setw(10) << event_count << '\r';
	    

    for (vector<McTrack>::const_iterator i = pevt->TrackList().begin(); i != pevt->TrackList().end(); ++i) {

      if (i->type() == 5) {

	// track parameters

	position  pos = getPosition(i->position());
	direction dir = getDirection(i->direction());
	double    t0  = i->t();

	const rotation R(dir);

	pos.rotate(R);
	
	for (vector<Hit>::const_iterator j = pevt->HitList().begin(); j != pevt->HitList().end(); ++j) {

	  if (i->id() == j->origin()) {
	    
	    try {
	      
	      const int    type   = j->type();
	      const double t1     = j->pure_dt();
	      const double npe    = j->pure_npe();
	      
	      const ars_address id = interface.getAddress(j->pm_id(),0);
	      
	      position  u  = interface.getPosition (id);
	      direction v  = interface.getDirection(id);
	    
	      u.rotate(R);
	      v.rotate(R);
	      
	      const double x   = u.x() - pos.x();
	      const double y   = u.y() - pos.y();
	      const double z   = u.z() - pos.z();
	      
	      const double r   = sqrt(x*x + y*y);

	      if (r >= RMIN && r <= RMAX) {
	      
		const double t1p = t0  +  (z  +  r * getTanThetaC()) * C_INVERSE;
		
		const double ct0 = 1.0 / getIndexOfRefraction();
		const double st0 = sqrt(1.0 - ct0*ct0);
		
		const double dx  = x * st0 / r;
		const double dy  = y * st0 / r;
		const double dz  = ct0;
		
		const double dot = dx*v.dx() + dy*v.dy() + dz*v.dz();
		
		if(type == +5 ||
		   type == -5 ||
		   type == +3 ||
		   type == -3) {
		  xmap[type]->Fill(t1 - t1p, npe);
		  //ymap[type]->Fill(dot,      npe);
		  ymap[type]->Fill(dot,      t1 - t1p);
		}
	      }
	    }
	    catch(const exception &error) {
	      cerr << error.what() << endl;
	      return 2;
	    }
	  }	
	}

	/*
	typedef TriggerInterfaceModel::alignment type;

	for (type::const_iterator j = interface.begin<type>(); j != interface.end<type>(); ++j) {

	  position  u  = (const position&)  j->second;
	  direction v  = (const direction&) j->second;
	    
	  u.rotate(R);
	  v.rotate(R);
	      
	  const double x   = u.x() - pos.x();
	  const double y   = u.y() - pos.y();
	  const double z   = u.z() - pos.z();
	  
	  const double r   = sqrt(x*x + y*y);
	  
	  if (r >= RMIN && r <= RMAX) {
		
	    const double ct0 = 1.0 / getIndexOfRefraction();
	    const double st0 = sqrt(1.0 - ct0*ct0);
	    
	    const double dx  = x * st0 / r;
	    const double dy  = y * st0 / r;
	    const double dz  = ct0;
	    
	    const double dot = dx*v.dx() + dy*v.dy() + dz*v.dz();

	    h0.Fill(dot, 1.0);
	  }
	}
	*/
      }
    }
  }
  cerr << endl; 

  /*
  for (map<int,TH1*>::iterator i = xmap.begin(); i != xmap.end(); ++i) {

    double z = i->second->GetSumOfWeights();

    if (z != 0.0) {

      for (int j = 1; j <= i->second->GetNbinsX(); ++j) {

	double y = i->second->GetBinContent(j);

	i->second->SetBinContent(j,y/z);
      }
    }
  }


  for (map<int,TH1*>::iterator i = ymap.begin(); i != ymap.end(); ++i) {

    for (int j = 1; j <= i->second->GetNbinsX(); ++j) {

      double y = i->second->GetBinContent(j);
      double z = h0.        GetBinContent(j);

      if (z != 0.0)
	i->second->SetBinContent(j,y/z);
      else
	i->second->SetBinContent(j,0.0);
    }
  }
  */

  {
    const double y = h0.GetSumOfWeights();

    h0.Scale(1.0/y);
  }

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