#ifndef JSHOWERENERGYPREFIT_INCLUDE
#define JSHOWERENERGYPREFIT_INCLUDE

#include <string>
#include <iostream>
#include <set>
#include <vector>
#include <algorithm>
#include <memory>
#include <math.h>

#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/JShowerEnergyRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JShowerNPEHit.hh"
#include "JFit/JShowerNPE.hh"

#include "JAAnet/JAAnetToolkit.hh"

#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerEnergyPrefitParameters_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::JSimplex;

  /**
   * class to handle the third step of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerEnergyPrefit : 
    public JShowerEnergyPrefitParameters_t,
    public JRegressor<JEnergy, JSimplex>
  {

    typedef JRegressor<JEnergy, JSimplex>  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
     */
    JShowerEnergyPrefit(const JShowerEnergyPrefitParameters_t&  parameters,
			const JModuleRouter&           router,
			const JSummaryRouter&          summary,
			const std::string              pdfFile,
			const int                      debug = 0):
      JShowerEnergyPrefitParameters_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->estimator.reset(getMEstimator(parameters.mestimator));
    }
    

    /**                                                                                                 
     * Auxiliary class for energy estimation.                                                
     */
    struct JResult {
      /**                                                                                                      
       * Constructor.                                                                      
       *                                                                                     
       * \param  x               Energy [log(E/GeV)]                                                          
       * \param  chi2            chi2                                                   
       */
      JResult(const JEnergy& x    = 0.0,
	      const double   chi2 = std::numeric_limits<double>::max())
      {
	this->x    = x;
	this->chi2 = chi2;
      }

      /**   
       * Type conversion.     
       * 
       * \return                 true if valid chi2; else false                                
       */
      operator bool() const
      {
	return chi2 != std::numeric_limits<double>::max();
      }

      JEnergy x;      //!< Energy                                                          
      double  chi2;   //!< chi2                                                                   
    };


    /**
     * Declaration of the member function that actually performs the reconstruction
     *
     * \param event which is a JDAQEvent
     * \param in input fits, which should contain a vertex fit
     */ 
    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;

	vector<JShowerNPEHit> data;
	vector<JShowerNPEHit> buffer;

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

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

	  if (match(*i)) {
	    top.insert(i->getPMTIdentifier());
	  }
	}

	JDetectorSubset_t subdetector(detector, sh.getPosition(), roadWidth_m);
	
	for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {

	  for (size_t i = 0; i != module->size(); ++i) {
	    
	    const size_t count  = top.count(JDAQPMTIdentifier(module->getID(), i));

	    const double rate_Hz = summary.getRate(JDAQPMTIdentifier(module->getID(), i));

	    data.push_back(JShowerNPEHit(this->getNPE(module->getPMT(i), rate_Hz), count));
	  }
	}

	const int NDF = data.size() - 1;

	if (NDF >= 0) {
	  
	  // 5-point search between given limits
	  const int      N  =  5;

	  JResult result[N];

	  for (int i = 0; i != N; ++i) {

	    result[i].x = log10(Emin_GeV + i * (Emax_GeV - Emin_GeV) / (N-1));

	  }
          
	  do{
  
	    int j = 0;

	    for (int i = 0; i != N; ++i) {
	      
	      if (!result[i]) {

		result[i].chi2 = (*this)(result[i].x, data.begin(), data.end());
	      
	      }
	      
	      if (result[i].chi2 < result[j].chi2) {
		j = i;
	      }
	    }	    

	    // squeeze range

           switch (j) {

           case 0:
             result[4] = result[1];
             result[2] = JResult(0.5 * (result[0].x + result[4].x));
             break;

           case 1:
             result[4] = result[2];
             result[2] = result[1];
             break;

           case 2:
             result[0] = result[1];
             result[4] = result[3];
             break;

           case 3:
             result[0] = result[2];
             result[2] = result[3];
             break;

           case 4:
             result[0] = result[3];
             result[2] = JResult(0.5 * (result[0].x + result[4].x));
             break;
           }
	   
           result[1] = JResult(0.5 * (result[0].x + result[2].x));
	   result[3] = JResult(0.5 * (result[2].x + result[4].x));
	   
	  } while (result[4].x - result[0].x > resolution);
      
	
	  if (result[1].chi2 != result[3].chi2) {
	    
	    result[2].x += 0.25 * (result[3].x - result[1].x) * (result[1].chi2 - result[3].chi2) / (result[1].chi2 + result[3].chi2 - 2*result[2].chi2);

	    result[2].chi2  = (*this)(result[2].x, data.begin(), data.end());
	  
	  }
         
	  //	  const double chi2 = result[2].chi2;  // this is not used because fits are sorted wrt the position fit
	  const double E    = result[2].x.getE();

	  JShower3E sh_fit(sh.getPosition(), sh.getDirection(), sh.getT(), E);
	
	  // the fits of this reco step are sorted wrt the previous reco step 
	  // because otherwise it tends to degrade the position reco performances
	  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERENERGYPREFIT), sh_fit, shower->getQ(), 
	  		       shower->getNDF(), sh_fit.getE())); 

	}
      }
      
      return out;

    }

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

#endif