#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 "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
 */

namespace JRECONSTRUCTION
{

  using JLANG::JSharedPointer;
  using JDETECTOR::JModuleRouter;

  /**
   * 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;

    /**
     * Parameterized constructor
     *
     * \param parameters    struct that holds default-optimized parameters for the reconstruction.
     * \param router        module router, this is built via detector file.
     * \param debug         debug
     */

    JShowerPrefit(const JShowerPrefitParameters_t&  parameters,
		  const JModuleRouter&              router,
		  const int                         debug = 0):
      JShowerPrefitParameters_t(parameters),
      router(router)
    {}
    
    /**
     * Declaration of the operator that performs the reconstruction
     *
     * \param event which is a JDAQEvent
     */

    JEvt operator()(const KM3NETDAQ::JDAQEvent &event) const
    {

      using namespace std;
      using namespace JPP;

      const double STANDARD_DEVIATIONS       = 3.0;

      typedef JEstimator<JPoint4D> JEstimator_t;

      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));
      const JMatch3G<hit_type> match3G(roadWidth_m, TMaxLocal_ns);  // causality relation for showers 
  
      //define empty vectors to be filled with JHitL0 and JHitL1
      buffer_type dataL0;
      buffer_type dataL1;
  
      JEvt out;   

      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));

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

	if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxLocal_ns)) == dataL1.end()) {
	  dataL0.push_back(*i); 
        } 
      }
      
      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(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;
	    }
	  }	
	}

	if (NDF >= 0 && chi2 > numeric_limits<double>::lowest()) {
        
	  const JGEOMETRY3D::JTrack3D fitResult(vx, JGEOMETRY3D::JVersor3Z());
	
	  out.push_back(getFit(JHistory(JSHOWERPREFIT), fitResult, getQuality(chi2, N), NDF));
	
	  out.rbegin()->setW(13, chi2);
	  out.rbegin()->setW(14, N);	       
	}
      }
   
      return out;
    
    }
        
    const JModuleRouter& router;

    /**
     * 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