#ifndef JSHOWERFIT_INCLUDE
#define JSHOWERFIT_INCLUDE

#include <string>
#include <iostream>
#include <set>
#include <vector>
#include <algorithm>
#include <memory>

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

#include "JTrigger/JHit.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.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 "JSupport/JSummaryRouter.hh"

#include "JFit/JShower3EZRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JGandalf.hh"

#include "JAAnet/JAAnetToolkit.hh"

#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerEnergyCorrection.hh"
#include "JReconstruction/JShowerFitParameters_t.hh"
#include "JReconstruction/JModuleL0.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"

#include "JDynamics/JCoverage.hh"

#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JGeometry3D/JShower3E.hh"

#include "JLang/JSharedPointer.hh"


/**
 * \author adomi
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {

  using KM3NETDAQ::JDAQEvent;
  using KM3NETDAQ::JDAQSummaryFrame;
  using JDETECTOR::JModuleRouter;
  using JDYNAMICS::coverage_type;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JShowerEnergyCorrection;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JFIT::JShower3EZ;
  using JFIT::JGandalf;
  /**
   * class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerFit : 
    public JShowerFitParameters_t,
    public JRegressor<JShower3EZ, JGandalf>
  {

    typedef JRegressor<JShower3EZ, JGandalf>  JRegressor_t;
    typedef JModuleL0                         module_type;
    typedef std::vector<module_type>          detector_type;
    using JRegressor_t::operator();
    
  public:


    /**
     * 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;
      detector_type data;
      coverage_type coverage;
    };

    /**
     * Parameterized constructor
     *
     * \param parameters   struct that holds default-optimized parameters for the reconstruction, available in $JPP_DATA. 
     * \param storage      storage
     * \param correct      energy correction
     * \param debug        debug
     */
    JShowerFit(const JShowerFitParameters_t&   parameters,
	       const storage_type&             storage,
	       const JShowerEnergyCorrection&  correct,
	       const int                       debug = 0):
      JShowerFitParameters_t(parameters),
      JRegressor_t(storage),
      correct(correct)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = debug;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
      JRegressor_t::Vmax_npe = parameters.Vmax_npe;
      JRegressor_t::MAXIMUM_ITERATIONS = 1000;
      JRegressor_t::EPSILON = 1e-3;
      JRegressor_t::EPSILON_ABSOLUTE = true;
      
      this->parameters.resize(3);

      this->parameters[0] = JShower3EZ::pDX();
      this->parameters[1] = JShower3EZ::pDY();
      this->parameters[2] = JShower3EZ::pE();

      this->estimator.reset(getMEstimator(parameters.mestimator));
    }

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

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

      const JBuildL0 <JHitR0>   buildL0;
      map<int, vector<JHitR0> > data;

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

	  buildL0(buffer, back_inserter(data[i->getModuleID()]));
        }
      }

      for (const auto& module : router.getReference()) {
	if (!module.empty()) {
	  input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID()), data[module.getID()]));
	}
      }

      return input;
    }

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

      JEvent event(JSHOWERDIRECTIONPREFIT);

      JEvt out;
      
      JEvt in = input.in;

      in.select(numberOfPrefits, qualitySorter);

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

      for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {

        vector<JPMTW0> data;

        const JPosition3D vertex(getPosition(*shower));
        const double      time = shower->getT();
        const double      distance = DMax_m  +  DStep_m * log10(shower->getE());
        const JRotation3D R(getDirection(*shower));
	
        for (const auto& module : input.data) {

          JPosition3D pos(module->getPosition());
          pos.sub(vertex);

          if(pos.getLength() < distance){

            for (size_t i = 0; i != module->size(); ++i) {
              if (module.getStatus(i)) {

		const double t1 = time + pos.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();

		struct {

                  bool operator()(const JHitR0& hit) const
                  {
                    return (hit.getPMT() == pmt && T_ns(hit.getT()));
                  }

                  const JTimeRange T_ns;
                  const size_t     pmt;

                } match = { JRegressor_t::T_ns + t1, i };

                JPMT pmt = module->getPMT(i);
		pmt.sub(vertex);
		pmt.rotate(R);
		
                data.push_back(JPMTW0(pmt, module.frame.getRate(i), count_if(module.begin(), module.end(), match)));
              }
            }
          }
        }
	    

	double chi2 = (*this)(JShower3EZ(JVertex3D(JVector3D(0,0,0), shower->getT()), JVersor3Z(), 
					 shower->getE()), data.begin(), data.end());
	  
	double NDF = getCount(data.begin(), data.end()) - this->parameters.size(); 

	JShower3E sh_fit(this->value.getPosition(), this->value.getDirection(),
			 this->value.getT(), correct(this->value.getE()));

	// check error matrix
	bool status = true;

	for (size_t i = 0; i != this->V.size(); ++i) {
	  if (std::isnan(this->V(i,i)) || this->V(i,i) < 0.0) {
	    status = false;
	  }
	}
	
	if (status) {
	  sh_fit.rotate_back(R);

	  sh_fit.add(vertex.getPosition());
	  
	  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT), sh_fit, getQuality(chi2), 
			       NDF, sh_fit.getE())); 
	  out.rbegin()->setV(this->V.size(), this->V);  
	  out.rbegin()->setW(JSHOWERFIT_ENERGY,  this->value.getE());  // Uncorrected Energy
	  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;
    }
    
    const JShowerEnergyCorrection&  correct;
  };
}

#endif