#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 "JAAnet/JAAnetToolkit.hh"

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

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

#include "JPhysics/KM3NeT.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, vcarretero
 */
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::JAbstractMinimiser;

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

    typedef JRegressor<JShower3EZ, JAbstractMinimiser>  JRegressor_t;
    using JRegressor_t::operator();
    vector<double> Ev;
    JPP::JOmega3D omega;

  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),
      omega(parameters.scanAngle_deg * JMATH::PI / 180.0),
      router(router),
      summary(summary)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = debug;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
      JRegressor_t::Vmax_npe = VMax_npe;
	  
      if (Emin_GeV > Emax_GeV || En <= 1) {
        THROW(JException, "Invalid energy input " << Emin_GeV << ' ' << Emax_GeV << ' ' << En);
      }

      const double base = std::pow((Emax_GeV / Emin_GeV), 1.0 / (En - 1));
      
      for (int i = 0; i != En; ++i) {
	      Ev.push_back(Emin_GeV * std::pow(base, i));
      }
    }

    /**
     * Declaration of the member function that actually performs the reconstruction
     *
     * \param event   DAQ event
     * \param in      input fits
     * \return        result 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) {

	const JVertex3D Vertex(getPosition(*shower),shower->getT());
	multiset<JDAQPMTIdentifier> top;
	
	const JFIT::JModel<JPoint4D> match(Vertex, DMax_m, JRegressor_t::T_ns);
	
	for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
	  if (match(*i)) {
	    top.insert(i->getPMTIdentifier());
	  }
	}

	const JDetectorSubset_t subdetector(detector, Vertex.getPosition(), DMax_m);
	
	vector<JPMTW0> data;
	for (JDetectorSubset_t::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
	  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());

	  JModule dom(*module);

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

	    if (getDAQStatus(frame, *module, i)            &&
		getPMTStatus(frame, *module, i)            &&
		frame[i].is_valid()                        &&
		!module->getPMT(i).has(PMT_DISABLE)) {

	      const JDAQPMTIdentifier id(module->getID(), i);

	      const double rate_Hz = summary.getRate(id);
	      const size_t count   = top.count(id);

	      data.push_back(JPMTW0(dom.getPMT(i), rate_Hz, count));
	    }
	  }
	}

	const JShower3EZ sh(JVertex3D(JVector3D(0,0,0), sh.getT()), JVersor3Z(), 1);

	for (JOmega3D_t::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
	  const JRotation3D R(*dir);
	  vector<double> chi2v(En, 0.0);
	  for (vector<JPMTW0>::const_iterator i = data.begin(); i != data.end(); ++i) {  
	    JPMTW0 pmt = *i;
	    pmt.rotate(R);				
	    JNPE_t::result_type H1 = (*this).getH1(sh, pmt);
	    JNPE_t::result_type H0 = (*this).getH0(pmt.getR());  // getH0 = Get background hypothesis value for time integrated PDF. 
	    const bool hit = pmt.getN() != 0;
	    for(size_t j=0;j!=chi2v.size();++j){
	      chi2v[j]+= getChi2(get_value(H0+Ev[j]*H1), hit); // H1 is linear with E
	    }
	  }
	  //Store only the lowest chi2 for a given direction
	  auto minChi2 = std::min_element(chi2v.begin(), chi2v.end());
	  JShower3E sh_fit(Vertex, JDirection3D(*dir), Ev[minChi2-chi2v.begin()]);
	  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERDIRECTIONPREFIT), sh_fit, getQuality(*minChi2), data.size(), sh_fit.getE()));
	}
      }
	
      return out;
    
    }
    const JModuleRouter&      router;
    const JSummaryRouter&     summary;
  };
}

#endif