#ifndef JSHOWERDIRECTIONPREFIT_INCLUDE
#define JSHOWERDIRECTIONPREFIT_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/JSimplex.hh"
#include "JFit/JGandalf.hh"

#include "JAAnet/JAAnetToolkit.hh"

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

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.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 JDETECTOR::JModuleRouter;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JFIT::JShower3EZ;
  using JFIT::JSimplex;
  using JFIT::JGandalf;

  /**
   * class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerDirectionPrefit : 
    public JShowerDirectionPrefitParameters_t,
    public JRegressor<JShower3EZ, JGandalf>
  {

    typedef JRegressor<JShower3EZ, JGandalf>  JRegressor_t;
    using JRegressor_t::operator();
    
  public:

    /**
     * Parameterized constructor
     *
     * \param parameters   struct that holds default-optimized parameters for the reconstruction, available in $JPP_DATA. 
     * \param router       module router, this is built via detector file. 
     * \param summary      summary router 
     * \param pdfFile      PDF file
     * \param debug        debug
     */
    JShowerDirectionPrefit(const JShowerDirectionPrefitParameters_t&   parameters,
			   const JModuleRouter&                        router,
			   const JSummaryRouter&                       summary,
			   const std::string                           pdfFile,
			   const int                                   debug = 0):
      JShowerDirectionPrefitParameters_t(parameters),
      JRegressor_t(pdfFile),
      router(router),
      summary(summary)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = 1;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
      JRegressor_t::Vmax_npe = VMax_npe;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;

      this->parameters.resize(2);
      
      this->parameters[0] = JShower3EZ::pDX();
      this->parameters[1] = JShower3EZ::pDY();
      
      this->estimator.reset(getMEstimator(parameters.mestimator));
    }

    /**
     * Declaration of the member function that actually performs the reconstruction
     *
     * \param event = JDAQEvent
     * \param in = input fits
     */ 
    JEvt operator()(const KM3NETDAQ::JDAQEvent& event, const JFIT::JEvt &in) 
    {
      using namespace std;
      using namespace JPP;

      typedef vector<JHitL0>  JDataL0_t;
      JBuildL0<JHitL0>        buildL0;

      JEvt out;

      const JDetector &detector = router.getReference();

      JDataL0_t dataL0;
      buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
  
      for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {

	JShower3E sh(getPosition(*shower), getDirection(*shower), shower->getT(), shower->getE());
	
	multiset<JDAQPMTIdentifier> top;

	double distance   = roadWidth_m;
	double max_angle  = MaxAngle_deg_lowE;
	double scan_angle = scanAngle_deg_lowE;
	
        if(sh.getE() > E_GeV) {
	  distance   = DMin_m;  // to speed up the code 
	  max_angle  = MaxAngle_deg_highE;
	  scan_angle = scanAngle_deg_highE;
	}

	const JFIT::JModel<JPoint4D> match(JPoint4D(sh.getPosition(), sh.getT()), distance, JRegressor_t::T_ns);

 	JVector3D start_dir(0,0,0);
	
	for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {

	  if (match(*i)) {

	    top.insert(i->getPMTIdentifier());

	    start_dir.add(JDirection3D(i->getPosition() - sh.getPosition()));	  

	  }
	}

	sh.setDirection(JDirection3D(start_dir));

	JDetectorSubset_t subdetector(detector, sh.getPosition(), distance);

	JOmega3D scan_directions(sh.getDirection(), 
				 JOmega3D::range_type(0, max_angle * PI/180.0), 
				 scan_angle * PI/180.0);

	for (JOmega3D_t::const_iterator dir = scan_directions.begin(); dir != scan_directions.end(); ++dir) {
	  
	  const JDirection3D conversion(*dir);	    
	  const JRotation3D R(conversion);
	
	  vector<JPMTW0> buffer;
	  
	  for (JDetectorSubset_t::iterator module = subdetector.begin();
	       module != subdetector.end(); ++module) {
	    
	    JModule dom(*module); 
	    dom.rotate(R);  

	    for (unsigned int i = 0; i != dom.size(); ++i) {

	      const JDAQPMTIdentifier id(dom.getID(), i);

	      const double rate_Hz = summary.getRate(id);

	      buffer.push_back(JPMTW0(dom.getPMT(i), rate_Hz, top.count(id)));
	    }      
	  }

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

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

	  sh_fit.add(sh.getPosition());
	  
	  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERDIRECTIONPREFIT), sh_fit, getQuality(chi2), 
			       NDF, sh_fit.getE())); 
	
	}
      }

      return out;

    }

    const JModuleRouter&      router;
    const JSummaryRouter&     summary;
  };
}

#endif