#ifndef  __JRECONSTRUCTION__JMUONSIMPLEX__
#define  __JRECONSTRUCTION__JMUONSIMPLEX__

#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>

#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"

#include "JTrigger/JHitR1.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"

#include "JFit/JFitToolkit.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JMEstimator.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonSimplexParameters_t.hh"

#include "JLang/JPredicate.hh"

#include "JDetector/JModuleRouter.hh"

#include "JDynamics/JCoverage.hh"

#include "JGeometry3D/JRotation3D.hh"

#include "Jeep/JMessage.hh"


/**
 * \author mdejong, gmaggi
 */

namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }

namespace JRECONSTRUCTION {

  using KM3NETDAQ::JDAQEvent;
  using JDETECTOR::JModuleRouter;
  using JFIT::JRegressor;
  using JFIT::JLine3Z;
  using JFIT::JSimplex;
  using JDYNAMICS::coverage_type;


  /**
   * Wrapper class to make intermediate fit of muon trajectory.
   *
   * The JMuonSimplex fit uses one or more start values (usually taken from the output of JMuonPrefit) and
   * produces new start values for subsequent fits (usually JMuonGandalf).\n
   * All hits of which the PMT position lies within a set road width (JMuonSimplexParameters_t::roadWidth_m) and
   * time is within a set window (JMuonSimplexParameters_t::TMin_ns, JMuonSimplexParameters_t::TMax_ns) around the Cherenkov hypothesis are taken.\n
   * In case there are multiple hits from the same PMT is the specified window,
   * the first hit is taken and the other hits are discarded.\n
   * The chi-squared is based on an M-estimator of the time residuals.\n
   */
  struct JMuonSimplex :
    public JMuonSimplexParameters_t,
    public JRegressor<JLine3Z, JSimplex>
  {
    typedef JRegressor<JLine3Z, JSimplex>        JRegressor_t;
    typedef JTRIGGER::JHitR1                     hit_type;
    typedef std::vector<hit_type>                buffer_type;      

    using JRegressor_t::operator();

    /**
     * Input data type.
     */
    struct input_type :
      public JDAQEventHeader
    {
      /**
       * Default constructor.
       */
      input_type()
      {}


      /**
       * Constructor.
       *
       * \param  header         header
       * \param  in             start values
       * \param  coverage       coverage
       */
      input_type(const JDAQEventHeader& header,
		 const JEvt&            in,
		 const coverage_type&   coverage) :
         JDAQEventHeader(header),
	 in(in),
	 coverage(coverage)
      {}

      JEvt          in;
      buffer_type   dataL0;
      buffer_type   dataL1;
      coverage_type coverage;
     };


    /**
     * Constructor 
     * 
     * \param  parameters    parameters
     * \param  debug         debug
     */
    JMuonSimplex(const JMuonSimplexParameters_t& parameters,
		 const int                       debug = 0) :
      JMuonSimplexParameters_t(parameters),
      JRegressor_t(parameters.sigma_ns)
    {
      using namespace JFIT;

      this->estimator.reset(new JMEstimatorLorentzian());

      JRegressor_t::debug              = debug;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;
    }
       

    /**
     * Get input data.
     *
     * \param  router        module router
     * \param  event         event
     * \param  in            start values
     * \param  coverage      coverage
     * \return               input data
     */
    input_type getInput(const JModuleRouter& router,
			const JDAQEvent&     event,
			const JEvt&          in,
			const coverage_type& coverage) const
    {
      using namespace std;
      using namespace JTRIGGER;
      using namespace KM3NETDAQ;

      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));

      input_type   input(event.getDAQEventHeader(), in, coverage);

      buffer_type& dataL0 = input.dataL0;
      buffer_type& dataL1 = input.dataL1;

      const JDAQTimeslice timeslice(event, true);
      JSuperFrame2D<JHit> buffer;

      for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {

        if (router.hasModule(i->getModuleID())) {

          buffer(*i, router.getModule(i->getModuleID()));

          if (useL0) {
            buildL0(buffer, back_inserter(dataL0));
          }
          {
            buildL2(buffer, back_inserter(dataL1));
          }
        }
      }

      return input;
    }


    /**
     * Fit function.
     *
     * \param  input         input data
     * \return               fit results
     */
    JEvt operator()(const input_type& input)
    {
      using namespace std;
      using namespace JFIT;
      using namespace JGEOMETRY3D;
      using namespace JLANG;

      JEvent event(JMUONSIMPLEX);

      JEvt out;

      const buffer_type& dataL0 = input.dataL0;
      const buffer_type& dataL1 = input.dataL1;

      // select start values

      JEvt in = input.in;

      in.select(numberOfPrefits, qualitySorter);

      if (!in.empty()) {
	in.select(JHistory::is_event(in.begin()->getHistory()));
      }

      buffer_type data;

      data.reserve(dataL0.size() + 
		   dataL1.size());
      
      for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {

	const JRotation3D     R (getDirection(*track));
	const JLine1Z         tz(getPosition (*track).rotate(R), track->getT());
	const JModel<JLine1Z> match(tz, roadWidth_m, JTimeRange(TMin_ns, TMax_ns));
	    
	data.clear();

	for (buffer_type::const_iterator i = dataL1.begin(); i != dataL1.end(); ++i) {

	  hit_type hit(*i);
	      
	  hit.rotate(R);

	  if (match(hit)) {
	    data.push_back(hit);
	  }
	}

	if (useL0) {

	  buffer_type::iterator __end = data.end();	

	  for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
	      
	    if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {

	      hit_type hit(*i);
	      
	      hit.rotate(R);
	      
	      if (match(hit)) {
		data.push_back(hit);
	      }
	    }
	  }
	}


	this->step.resize(5);
	
	this->step[0] = JLine3Z(JLine1Z(JVector3D(0.5, 0.0, 0.0), 0.0));
	this->step[1] = JLine3Z(JLine1Z(JVector3D(0.0, 0.5, 0.0), 0.0));
	this->step[2] = JLine3Z(JLine1Z(JVector3D(0.0, 0.0, 0.0), 1.0));
	this->step[3] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.005, 0.0));
	this->step[4] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.0, 0.005));
	
	const int NDF = getCount(data.begin(), data.end()) - this->step.size();

	if (NDF > 0) {

	  const double chi2 = (*this)(JLine3Z(tz), data.begin(), data.end());
	      
	  JTrack3D tb(this->value);
	      
	  tb.rotate_back(R);

	  out.push_back(getFit(JHistory(track->getHistory(), event()), tb, getQuality(chi2, NDF), NDF)); 

	  // set additional values

	  out.rbegin()->setW(track->getW());
          out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
          out.rbegin()->setW(JPP_COVERAGE_POSITION,    input.coverage.position);
  	}
      }

      // apply default sorter

      sort(out.begin(), out.end(), qualitySorter);

      copy(input.in.begin(), input.in.end(), back_inserter(out));

      return out;
    }
  };
}

#endif