#ifndef JSHOWERPREFIT_INCLUDE
#define JSHOWERPREFIT_INCLUDE

#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>

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

#include "JDynamics/JCoverage.hh"

#include "JTrigger/JHit.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch3G.hh"
#include "JTrigger/JBind2nd.hh"

#include "JTools/JPermutation.hh"

#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4DEstimator.hh"
#include "JFit/JEstimator.hh"
#include "JFit/JPoint4D.hh"

#include "JDetector/JModuleRouter.hh"

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

#include "JLang/JComparator.hh"

#include "JGeometry3D/JTrack3D.hh"

/**
 * \author adomi, gmaggi, vcarretero
 */

namespace JRECONSTRUCTION
{
  using KM3NETDAQ::JDAQEvent;
  using JLANG::JSharedPointer;
  using JDETECTOR::JModuleRouter;
  using JDYNAMICS::coverage_type;
  
  /**
   * class to handle first step of the shower reconstruction in ORCA:
   * it reconstructs the shower vertex, intended as the shower brightest point, as the barycenter of the hits
   */

  class JShowerPrefit : public JShowerPrefitParameters_t
  {
    
  public:
    
    typedef JTRIGGER::JHitR1                     hit_type;
    typedef std::vector<hit_type>                buffer_type;

    /**
     * Input data type.
     **/
    struct input_type :
      public JDAQEventHeader
    {
      /**
       * Default constructor.
       */
      input_type() 
      {}
      /**
       * Constructor.
       *
       * \param  header         header
       * \param  coverage      coverage
       */
      input_type(const JDAQEventHeader& header, const coverage_type& coverage) :
	JDAQEventHeader(header),
	coverage       (coverage)
      {}

      buffer_type  dataL0;
      buffer_type  dataL1;
      coverage_type coverage;
    };

    /**
     * Parameterized constructor
     *
     * \param parameters    struct that holds default-optimized parameters for the reconstruction.
     * \param debug         debug
     */

    JShowerPrefit(const JShowerPrefitParameters_t&  parameters,
		  const int                         debug = 0):
      JShowerPrefitParameters_t(parameters)
    {}
    
    /**
     * Get input data.
     *
     * \param  router        module router
     * \param  event         event
     * \param  coverage      coverage
     * \return               input data
     */
    input_type getInput(const JModuleRouter& router, const KM3NETDAQ::JDAQEvent& event, 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(), coverage);
      
      buffer_type& dataL0 = input.dataL0;
      buffer_type& dataL1 = input.dataL1;
  
      buffer_type buffer;

      
      buildL0(JDAQTimeslice(event, true),  router, back_inserter(buffer));  // true  => snapshot
      buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1));  // false => triggered
      
      copy(dataL1.begin(), dataL1.end(), back_inserter(dataL0));

      if(dataL1.size() <= numberOfL1){
	for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
	  if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxExtra_ns)) == dataL1.end()) {
	    dataL0.push_back(*i); 
	  } 
	}
      }
      return input;
    }
    /**
     * Fit function.
     *
     * \param  input         input data
     * \return               fit results
     */
    JEvt operator()(const input_type& input)
    {
      using namespace std;
      using namespace JPP;

      const double STANDARD_DEVIATIONS       = 3.0;

      typedef JEstimator<JPoint4D> JEstimator_t;
      
      JEvent event(JSHOWERPREFIT);
      
      JEvt out; 

      const buffer_type& dataL0 = input.dataL0;
      const buffer_type& dataL1 = input.dataL1;
      
      const JMatch3G<hit_type> match3G(DMax_m, TMaxExtra_ns);  // causality relation for showers

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

	buffer_type data(1, *root);
      
	JBinder2nd<hit_type> matching = JBind2nd(match3G, *root);
      
	for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
	
	  if(( root->getModuleIdentifier() != i->getModuleIdentifier() ) && matching(*i)){
	    data.push_back(*i);	    
	  }
	}
	
	buffer_type::iterator __end1 = clusterizeWeight(data.begin() + 1, data.end(), match3G);
      
	// 4D fit
	JEstimator_t fit;
	JPoint4D  vx;
	double    chi2 = numeric_limits<double>::max();
	int       NDF  = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS;
	int	  N    = getCount(data.begin(), __end1);

	if(NDF > 0){
	  if(distance(data.begin(), __end1) <= factoryLimit){
	  
	    double ymin = numeric_limits<double>::max();
      
	    buffer_type::iterator __end2 = __end1;

	    for (int n = 0; n <= numberOfOutliers && distance(data.begin(), __end2) > 
		   JEstimator_t::NUMBER_OF_PARAMETERS; ++n, --__end2) {
	
	      sort(data.begin() + 1, __end1, hit_type::compare);

	      do {
		try {
	    
		  fit(data.begin(), __end2);
		  
		  double y = getChi2(fit, data.begin(), __end2, sigma_ns);

		  if (y < ymin) {
		    ymin = y;
		    vx   = fit;
		    chi2 = ymin;
		    NDF  = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
		    N    = getCount(data.begin(), __end2);
		  }
		}
		catch(JException& error) { }
	  
	      } while (next_permutation(data.begin() + 1, __end2, __end1, hit_type::compare));
	
	      ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
	    }
	    
	  } else {

	    const int number_of_outliers = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS - 1;

	    buffer_type::iterator __end2 = __end1;

	    for (int n = 0; n <= number_of_outliers; ++n) {	  

	      try{	

		fit(data.begin(), __end2);
		vx   = fit;
		chi2 = getChi2(fit, data.begin(), __end2, sigma_ns);
		NDF  = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
		N    = getCount(data.begin(), __end2);
		
	      }
	      catch(const JException& error){ }

	      double                ymax = 0;
	      buffer_type::iterator imax = __end2;

	      for (buffer_type::iterator i = data.begin() + 1; i != __end2; ++i) {

		double y = getChi2(fit, *i, sigma_ns);

		if (y > ymax) {
		  ymax = y;
		  imax = i;
		}
	      }

	      if (ymax > STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {	
		--__end2;
		swap(*imax, *__end2);	  
	      } else {
		break;
	      }
	    }	
	  }
  
	  out.push_back(getFit(JHistory(JSHOWERPREFIT), JTrack3D(vx, JGEOMETRY3D::JVersor3Z()), getQuality(chi2, N), NDF));

	  // set additional values

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

	}
      }

      out.select(numberOfGrids, qualitySorter);

      size_t solutions = out.size();

      for(size_t i=0; i < solutions; i++){
	for(int x = -pos_grid_m; x < pos_grid_m + pos_step_m/2.; x += pos_step_m){
	  for(int y = -pos_grid_m; y < pos_grid_m + pos_step_m/2.; y += pos_step_m){
	    for(int z = -pos_grid_m; z < pos_grid_m + pos_step_m/2.; z += pos_step_m){
	      for(int t = -time_grid_ns; t < time_grid_ns + time_step_ns/2.; t += time_step_ns){
		if (x != 0 || y != 0 || z != 0 || t != 0) {
		  
		  out.push_back(getFit(event(),
				       JTrack3D(JPoint4D(getPosition(out[i]) + JPosition3D(x,y,z), out[i].getT()+t), JVersor3Z()),
				       0,
				       0));
		  
		  // set additional values

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

		}
	      }
	    }
	  }
	}
      }
      //sorting
      
      if (numberOfPrefits > 0) {
 
	// apply default sorter
      
	JFIT::JEvt::iterator __end = out.end();

	if (numberOfPrefits < out.size()) {

	  advance(__end = out.begin(), numberOfPrefits);

	  partial_sort(out.begin(), __end, out.end(), qualitySorter);

	  out.erase(__end, out.end());

	} else {

	  sort(out.begin(), __end, qualitySorter);
	}

      } else {
      
	sort(out.begin(), out.end(), qualitySorter);
      }

      return out;
    }
    
    /**
     * Auxiliary class to match hit to root hit.
     */
    struct match_t {
      /**
       * Constructor.
       *
       * \param  root        root hit
       * \param  TMax_ns     maximal time [ns]
       */
      match_t(const hit_type& root, const double TMax_ns) :
	root(root),
	TMax_ns(TMax_ns)
      {}

      /**
       * Test match.
       *
       * \param  hit         hit
       * \return             true if mach; else false
       */
      bool operator()(const hit_type& hit) const
      {
	return root.getModuleID() == hit.getModuleID() && fabs(root.getT() - hit.getT()) <= TMax_ns;
      }

      const hit_type& root;
      double TMax_ns;
    };
  };
}

#endif