#ifndef __JRECONSTRUCTION__JMUONENERGY__
#define __JRECONSTRUCTION__JMUONENERGY__

#include <string>
#include <istream>
#include <ostream>
#include <iomanip>
#include <set>
#include <vector>
#include <algorithm>
#include <cmath>

#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"

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

#include "JDynamics/JCoverage.hh"

#include "JTrigger/JHitR0.hh"
#include "JTrigger/JBuildL0.hh"

#include "JSupport/JSummaryRouter.hh"

#include "JFit/JLine1Z.hh"
#include "JFit/JEnergyRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JModel.hh"
#include "JFit/JNPEHit.hh"
#include "JFit/JEnergy.hh"

#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonEnergyParameters_t.hh"
#include "JReconstruction/JEnergyCorrection.hh"
#include "JReconstruction/JModuleL0.hh"

#include "JPhysics/JGeane.hh"


/**
 * \author mdejong, azegarelli, scelli
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {

  using KM3NETDAQ::JDAQEvent;
  using KM3NETDAQ::JDAQSummaryFrame;
  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JModule;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JDYNAMICS::coverage_type;


  /**
   * Auxiliary class to to determine muon energy.
   */
  class JMuonEnergy : 
    public JMuonEnergyParameters_t,
    public JRegressor<JEnergy>
  {
  public:
    typedef JRegressor<JEnergy>                  JRegressor_t;
    typedef JModuleL0                            module_type;
    typedef std::vector<module_type>             detector_type;

    using JRegressor_t::operator();

    /**
     * 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)
      {}

      JEvt          in;
      detector_type data;
      coverage_type coverage;
    };


    /**
     * Constructor.
     *
     * \param  parameters    parameters
     * \param  storage       storage
     * \param  correct       energy correction
     * \param  debug         debug
     */
    JMuonEnergy(const JMuonEnergyParameters_t& parameters,
		const storage_type&            storage,
		const JEnergyCorrection&       correct,
		const int                      debug = 0):
      JMuonEnergyParameters_t(parameters),
      JRegressor_t(storage),
      correct     (correct)
    {
      using namespace JPP;
      
      if (this->getRmax() < roadWidth_m) {
        roadWidth_m = this->getRmax();
      }

      JRegressor_t::debug  = debug;
      JRegressor_t::T_ns.setRange(TMin_ns, TMax_ns);

      this->estimator.reset(getMEstimator(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) const
    {
      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;
    }
    

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


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

      JEvent event(JMUONENERGY);

      JEvt out;

      // select start values

      JEvt in  = input.in;

      in.select(numberOfPrefits, qualitySorter);

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

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

	const JRotation3D  R (getDirection(*track));
	const JLine1Z      tz(getPosition (*track).rotate(R), track->getT());

        double zmin = numeric_limits<double>::lowest();

        if (track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0) {
          zmin = this->ZMin_m;
        }

	vector<JNPEHit> data;

        for (const auto& module : input.data) {

          JPosition3D pos(module->getPosition());

          pos.transform(R, tz.getPosition());

          if (pos.getX() <= roadWidth_m) {

            const double z1 = pos.getZ()  -  pos.getX() / getTanThetaC();
	    const double t1 = tz .getT()  + (pos.getZ() + pos.getX() * getKappaC()) * getInverseSpeedOfLight();

	    if (z1 >= zmin) {

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

		if (module.getStatus(i)) {

		  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.transform(R, tz.getPosition());

		  data.push_back(JNPEHit(this->getNPE(pmt, module.frame.getRate(i)), count_if(module.begin(), module.end(), match)));
		}
	      }
	    }
	  }
	}

	const int NDF = distance(data.begin(), data.end()) - 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 = EMin_log + i * (EMax_log - EMin_log) / (N-1);
	  }
	  
	  map<double, double> buffer;
	  
	  do {
	    
	    int j = 0;
	    
	    for (int i = 0; i != N; ++i) {
	      
	      if (!result[i]) {
		
		const JEnergy x    = result[i].x;
		const double  chi2 = (*this)(x, data.begin(), data.end());
		
		result[i].chi2 = chi2;
		buffer[chi2]   = x.getE();
	      }
	      
	      if (result[i].chi2 < result[j].chi2) {
		j = i;
	      }
	    }  
	    
	    
	    for (int i = 0; i != N; ++i) {
	      DEBUG(' ' << FIXED(5,2) << result[i].x << ' ' << FIXED(9,3) << result[i].chi2);
	    }
	    DEBUG(endl);
	    
	    // squeeze range
	    
	    switch (j) {
	      
	    case 0:
	      result[4] = result[1];
	      result[2] = result[0];
	      result[0] = JResult(2*result[2].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] = result[4];
	      result[4] = JResult(2*result[2].x - result[0].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;
	  const double E    = result[2].x.getE();
	  
	  // calculate additional variables

	  double Emin = numeric_limits<double>::max();
	  double Emax = numeric_limits<double>::lowest();
	  
	  for (map<double, double>::const_iterator i = buffer.begin(); i != buffer.end() && i->first <= chi2 + 1.0; ++i) {
	    if (i->second < Emin) { Emin = i->second; }
	    if (i->second > Emax) { Emax = i->second; }
	  }
	  
	  const double mu_range   = gWater(E);           // range of a muon with the reconstructed energy
	  
	  double noise_likelihood = 0.0;                 // log-likelihood of every hit being from K40   
	  int    number_of_hits   = 0;                   // number of hits selected for JEnergy
	  
	  for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
	    noise_likelihood += log10(getP(i->getY0(), i->getN()));       // probability of getting the observed multiplicity with K40 
	    number_of_hits   += i->getN();                                // hit multiplicity
	  }

	  JFit fit = *track;

	  fit.push_back(event());

	  // set corrected energy

	  fit.setE(correct(E));

 	  out.push_back(fit);

	  // set additional values

	  out.rbegin()->setW(track->getW());
	  out.rbegin()->setW(JENERGY_ENERGY,            E);
	  out.rbegin()->setW(JENERGY_CHI2,              chi2);
	  out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, mu_range);
	  out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD,  noise_likelihood);
	  out.rbegin()->setW(JENERGY_NDF,               NDF);
	  out.rbegin()->setW(JENERGY_NUMBER_OF_HITS,    number_of_hits);
	  out.rbegin()->setW(JENERGY_MINIMAL_ENERGY,    Emin);
	  out.rbegin()->setW(JENERGY_MAXIMAL_ENERGY,    Emax);
          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;
    }

    JEnergyCorrection correct;
  };
}

#endif