#include <string>
#include <iostream>
#include <strstream>
#include <stdexcept>
#include <fstream>
#include <iomanip>
#include <map>
#include <algorithm>

#include "TFile.h"
#include "TH1.h"

#include "jeep/parser.hh"
#include "jeep/properties.hh"
#include "antcc/Event.hh"
#include "trigger/cluster.hh"
#include "trigger/montecarlo_helper.hh"
#include "nr/prob.hh"
#include "nr/random.hh"

#include "prefit.hh"
#include "monopole.hh"

#include "TriggeredEvent.hh"
#include "RHit.hh"
#include "REvent.hh"


int main(int argc, char **argv)
{
  using namespace std;
  using namespace TRIGGER;
  using namespace TRIGGER_IO;
  using namespace CLUSTER;
  using namespace PREFIT;
  using namespace MONOPOLE;

  string   inputFile;
  string   inputHistos;
  //string   muonHistos;
  string   detectorFile;
  size_t   numberOfEvents;
  double   grid;
  unsigned numberOfLCMs;
  float    maxBeta;
  float    minBeta;
  float    maxRange;
  float    maxRoadWidth;
  string   outputFile;
  bool     debug;

  Parser<> zap;

  zap['f'] = make_field(inputFile);
  zap['H'] = make_field(inputHistos);
  //zap['u'] = make_field(muonHistos);
  zap['a'] = make_field(detectorFile)   = "d10_c00_s00.det";
  zap['n'] = make_field(numberOfEvents) = 1;
  zap['g'] = make_field(grid)           = 10.0;  // [deg]
  zap['P'] = make_field(numberOfLCMs)   = 7; // BvR
  zap['b'] = make_field(maxBeta)        = 1;
  zap['B'] = make_field(minBeta)        = 1;
  zap['z'] = make_field(maxRange)       = 24;
  zap['R'] = make_field(maxRoadWidth)   = 90;
  zap['o'] = make_field(outputFile)     = "fit.evt";
  zap['d'] = make_field(debug);

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

  // initialise minuit
  mninit(debug);

  MonteCarloHelper interface(detectorFile.c_str());

  // set minimum correlation speed (for slowly moving particles): BvR
  setMaxBeta(maxBeta);
  setMinBeta(minBeta);
  setMaxRangeTime(maxRange);
  setMaxTransverseDistance(maxRoadWidth);

  cout << "Grid                        : " << grid << endl
       << "Minimum number of LCMs      : " << numberOfLCMs << endl
       << "Maximum Beta inversed       : " << getMaxBetaInverse() << endl
       << "Minimum Beta inversed       : " << getMinBetaInverse() << endl
       << "Max. range time             : " << getMaxRangeTime() << endl
       << "Min. range time             : " << getMinRangeTime() << endl
       << "Max. transverse distance    : " << getMaxTransverseDistance() << endl;

  interface.center(inputFile.c_str());

  ifstream fin(inputFile.c_str());
  ofstream fout(outputFile.c_str(), ios::out);
  TFile histos(inputHistos.c_str());
  //TFile muhistos(muonHistos.c_str());

  for (string buffer; getline(fin,buffer) && fout << buffer << endl && buffer != "end_event:"; ) {}

  //std::vector<direction> omega;
  Omega omega(grid);

  REvent event;

  TH1D *dn55 = (TH1D*) histos.Get("dn55");
  TH1D *dn60 = (TH1D*) histos.Get("dn60");
  TH1D *dn65 = (TH1D*) histos.Get("dn65");
  TH1D *dn70 = (TH1D*) histos.Get("dn70");
  TH1D *dn75 = (TH1D*) histos.Get("dn75");
  //TH1D *dnmu = (TH1D*) muhistos.Get("aresraw");
  //dnmu->Smooth(200);

  for (size_t eventCounter = 0; eventCounter < numberOfEvents && fin; ++eventCounter) { // BvR event loop

    event.Streamer(fin);
    if(!fin)  break;

    cerr << "event: " << setw(8) << event.eventNumber() << '\r';
    cout << "event: " << setw(8) << event.eventNumber() << endl;

    // event selection

    if (event.trigger() == 0)
      {
	//cout << "No trigger " << endl;
	event.Streamer(fout);	
	continue;
      }
    
    vector<RHit> data;
    
    for(vector<RawHit>::const_iterator i = event.triggeredRawHitList().begin(); i != event.triggeredRawHitList().end(); ++i)
      {
	const header& id = interface.getHeader(i->pm_id());
	
	data.push_back(RHit(id,
			    interface.getPosition(id),
			    *i,
			    25.0));
      }
    
    // BvR: count number of different LCMs
    set<unsigned> lcms;
    lcms.clear();
    for(vector<RHit>::const_iterator i = data.begin(); i != data.end(); i++)
      {
	lcms.insert(i->lcm_id());
      }
    if (lcms.size() < numberOfLCMs)
      {
	event.Streamer(fout);
	continue;
      }
    
    //double fmin = DBL_MAX;

    //omega.push_back(Dir);

    for (Omega::const_iterator i_ = omega.begin(); i_ != omega.end(); ++i_) { //BvR loop dir muon

      const direction& dir = *i_;
      const direction  dir0(0,0);
      
      rotation R(dir);
      
      for (vector<RHit>::iterator i = data.begin(); i != data.end(); ++i)
        i->rotate(R);
      
      
      try {
	
	//----------
	// MUON FIT
	//----------
	
	double BETA = 1;
	//hist = dnmu;
	//sigma = 0.261;
	sigma = 0.478;
	//cout << BETA << ' '  << hist->GetBinContent(2000) << ' ' << sigma << endl;
	
	// prefit
	
	prefit1Zmu pfitmu(data.begin(),data.end());
	    
	int n = (int) data.size();
	/*
	event.prefitmu.insert(make_pair(n,REvent::track1d(pfitmu.x(),
							  pfitmu.y(),
							  pfitmu.z(),
							  dir.theta(),
							  dir.phi(),
							  pfitmu.t(),
							  BETA,
							  pfitmu.chi2ndf_)));
	*/
	fitted_vtrack ft1mu = minuitmu(vtrack(pfitmu,dir0,BETA),
				       data.begin(), 
				       data.end(), 
				       "CALL FCN",
				       "FIX 6",
				       "MINIMIZE 10000",
				       NULL);
	
	if( (180/PI) * acos( dot( direction(ft1mu.theta(),ft1mu.phi()), dir0 ) ) <= 1.4 * grid )
	  { // if final fit
	    
	    ft1mu.rotate_back(R);
	    
	    event.fitmu.insert(make_pair(n,REvent::track3d(ft1mu.a(),
							   ft1mu.b(),
							   ft1mu.theta(),
							   ft1mu.phi(),
							   ft1mu.t0(),
							   ft1mu.beta(),
							   ft1mu.fmin(),
							   0)));
	  } // if final fit muon
	
	
	//#ifdef NONONONONONO	 
	//----------
	// MONOP FIT
	//----------
	
	
	// prefit
	
	prefit1Z pfit(data.begin(),data.end());
	
	n = (int) data.size();
	/*
	event.prefit.insert(make_pair(n,REvent::track1d(pfit.x(),
							pfit.y(),
							pfit.z(),
							dir.theta(),
							dir.phi(),
							pfit.t(),
							pfit.beta_,
							pfit.chi2ndf_)));
	*/
	BETA  = 0.55;
	hist  = dn55;
	//sigma = 0.1758;
	sigma = 0.394;
	//cout << BETA << ' '  << hist->GetBinContent(2000) << ' ' << sigma << endl;
	    
	fitted_vtrack ft1 = minuit(vtrack(pfit,dir0,BETA),
				   data.begin(), 
				   data.end(), 
				   "CALL FCN",
				   "FIX 6",
				   "MINIMIZE 10000",
				   NULL);
	
	if( (180/PI) * acos( dot( direction(ft1.theta(),ft1.phi()), dir0 ) ) <= 1.4 * grid )
	  { // if final fit
	    
	    ft1.rotate_back(R);
	    
	    event.fit.insert(make_pair(n,REvent::track3d(ft1.a(),
							 ft1.b(),
							 ft1.theta(),
							 ft1.phi(),
							 ft1.t0(),
							 ft1.beta(),
							 ft1.fmin(),
							 0)));
	  } // if final fit monop
	
	BETA  = 0.60;
	hist  = dn60;
	//sigma = 0.3069;
	sigma = 0.693;
	//cout << BETA << ' '  << hist->GetBinContent(2000) << ' ' << sigma << endl;
	    
	ft1 = minuit(vtrack(pfit,dir0,BETA),
		     data.begin(), 
		     data.end(), 
		     "CALL FCN",
		     "FIX 6",
		     "MINIMIZE 10000",
		     NULL);
	
	if( (180/PI) * acos( dot( direction(ft1.theta(),ft1.phi()), dir0 ) ) <= 1.4 * grid )
	  { // if final fit
	    
	    ft1.rotate_back(R);
	    
	    event.fit.insert(make_pair(n,REvent::track3d(ft1.a(),
							 ft1.b(),
							 ft1.theta(),
							 ft1.phi(),
							 ft1.t0(),
							 ft1.beta(),
							 ft1.fmin(),
							 0)));
	  } // if final fit monop
	
	BETA  = 0.65;
	hist  = dn65;
	//sigma = 0.4274;
	sigma = 0.978;
	//cout << BETA << ' '  << hist->GetBinContent(2000) << ' ' << sigma << endl;
	
	ft1 = minuit(vtrack(pfit,dir0,BETA),
		     data.begin(), 
		     data.end(), 
		     "CALL FCN",
		     "FIX 6",
		     "MINIMIZE 10000",
		     NULL);
	
	if( (180/PI) * acos( dot( direction(ft1.theta(),ft1.phi()), dir0 ) ) <= 1.4 * grid )
	  { // if final fit
	    
	    ft1.rotate_back(R);
	    
	    event.fit.insert(make_pair(n,REvent::track3d(ft1.a(),
							 ft1.b(),
							 ft1.theta(),
							 ft1.phi(),
							 ft1.t0(),
							 ft1.beta(),
							 ft1.fmin(),
							 0)));
	  } // if final fit monop
	
	BETA  = 0.70;
	hist  = dn70;
	//sigma = 0.5320;
	sigma = 1.25;
	//cout << BETA << ' '  << hist->GetBinContent(2000) << ' ' << sigma << endl;
	    
	ft1 = minuit(vtrack(pfit,dir0,BETA),
		     data.begin(), 
		     data.end(), 
		     "CALL FCN",
		     "FIX 6",
		     "MINIMIZE 10000",
		     NULL);
	
	if( (180/PI) * acos( dot( direction(ft1.theta(),ft1.phi()), dir0 ) ) <= 1.4 * grid )
	  { // if final fit
	    
	    ft1.rotate_back(R);
	    
	    event.fit.insert(make_pair(n,REvent::track3d(ft1.a(),
							 ft1.b(),
							 ft1.theta(),
							 ft1.phi(),
							 ft1.t0(),
							 ft1.beta(),
							 ft1.fmin(),
							 0)));
	  } // if final fit monop
	
	BETA  = 0.75;
	hist  = dn75;
	//sigma = 0.6700;
	sigma = 1.32;
	//cout << BETA << ' '  << hist->GetBinContent(2000) << ' ' << sigma << endl;
	    
	ft1 = minuit(vtrack(pfit,dir0,BETA),
		     data.begin(), 
		     data.end(), 
		     "CALL FCN",
		     "FIX 6",
		     "MINIMIZE 10000",
		     NULL);
	
	if( (180/PI) * acos( dot( direction(ft1.theta(),ft1.phi()), dir0 ) ) <= 1.4 * grid )
	  { // if final fit
	    
	    ft1.rotate_back(R);
	    
	    event.fit.insert(make_pair(n,REvent::track3d(ft1.a(),
							 ft1.b(),
							 ft1.theta(),
							 ft1.phi(),
							 ft1.t0(),
							 ft1.beta(),
							 ft1.fmin(),
							 0)));
	  } // if final fit monop
	
	//#endif	
	
      } // if try
      catch(exception& error) {
	cerr << error.what() << endl;
      }
      
      for (vector<RHit>::iterator i = data.begin(); i != data.end(); ++i)
	i->rotate_back(R);
      
    } // BvR loop dir muon
    
    //cout << "FMIN " << fmin << ' ' << data.size() << ' ' << fmin/(2*data.size()) << endl;
    
    //omega.pop_back();
    
    event.Streamer(fout);
  } // BvR event loop
  cerr << endl;
  
  fout.close();
}