#ifndef __JRECONSTRUCTION__JMUONSTART__
#define __JRECONSTRUCTION__JMUONSTART__

#include <string>
#include <istream>
#include <ostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#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 "JTrigger/JHitL0.hh"
#include "JTrigger/JHitR1.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 "JFit/JMEstimator.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonStartParameters_t.hh"
#include "JReconstruction/JStart.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JK40Rates.hh"

#include "JLang/JComparator.hh"
#include "JLang/JPredicate.hh"


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

  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JModule;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JPHYSICS::JK40Rates;
  using JPHYSICS::JRateL1_t;
  using KM3NETDAQ::JDAQSummaryFrame;

  /**
   * Get probability of given response in optical module due to random background.
   *
   * \param  module               module
   * \param  frame                summary frame
   * \param  R_Hz                 multiples rates [Hz]
   * \param  T_ns                 time window [ns]
   * \param  top                  list of indentifiers of PMTs with hit in given module
   * \return                      probability
   */
  inline double getProbability(const JModule&            module,
                               const JDAQSummaryFrame&   frame,
                               const JRateL1_t&          R_Hz,
                               const double              T_ns,
                               const std::multiset<int>& top)
  {
    using namespace std;
    using namespace JPP;
    using namespace KM3NETDAQ;

    size_t N = 0;         // number of active PMTs
    size_t M = 0;         // multiplicity
    double R = 0.0;       // total rate

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

      if (getDAQStatus(frame, module, i)            &&
          getPMTStatus(frame, module, i)            &&
          !module.getPMT(i).has(PMT_DISABLE)) {

        N += 1;
        M += top.count(i);
        R += frame.getRate(i);
      }
    }

    if (N > 0)
      return JFIT::getProbability(N, M, JK40Rates(R/N, R_Hz), T_ns);
    else
      return (M == 0 ? 1.0 : 0.0);
  }


  /**
   * Auxiliary class to determine start and end position of muon trajectory.
   */
  class JMuonStart : 
    public JMuonStartParameters_t,
    public JRegressor<JEnergy>
  {
  public:
    typedef JRegressor<JEnergy>                  JRegressor_t;
    typedef JTRIGGER::JHitL0                     hit_type;
    typedef std::vector<hit_type>                buffer_type;

    using JRegressor_t::operator();
    
    /**
     * Constructor.
     *
     * \param  parameters   parameters
     * \param  router       module  router
     * \param  summary      summary router
     * \param  pdfFile      PDF file
     * \param  rates_Hz     K40 rates [Hz]
     * \param  debug        debug
     */
    JMuonStart(const JMuonStartParameters_t& parameters,
	       const JModuleRouter&          router,
	       const JSummaryRouter&         summary,
	       const std::string&            pdfFile,
	       const JK40Rates&              rates_Hz,
	       const int                     debug = 0):
      JMuonStartParameters_t(parameters),
      JRegressor_t(pdfFile),
      router(router),
      summary(summary),
      rates_Hz(rates_Hz),
      debug(debug)
    {
      JRegressor_t::debug = debug;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);

      if (this->getRmax() < parameters.roadWidth_m) {
        roadWidth_m = this->getRmax();
      }
    }


    /**
     * Fit function.
     *
     * \param  event         event
     * \param  in            start values
     * \return               fit results
     */
    JEvt operator()(const KM3NETDAQ::JDAQEvent& event, const JEvt& in)
    {
      using namespace std;
      using namespace JPP;

      const JBuildL0<hit_type>  buildL0;

      buffer_type dataL0;

      buildL0(event, router, true, back_inserter(dataL0));

      return (*this)(dataL0, in);
    } 


    /**
     * Fit function.
     *
     * \param  data          dataL0
     * \param  in            start values
     * \return               fit results
     */
    JEvt operator()(const buffer_type& data, const JEvt& in)
    {
      using namespace std;
      using namespace JPP;

      JEvt out;

      const JDetector& detector = router.getReference();

      const JStart start(Pmin1, Pmin2, Nmax2);  
      
      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());
        const JFIT::JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
	
	map<int, multiset<int> > top;
	
	for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
	  
	  hit_type hit(*i);
	  
	  hit.rotate(R);
	  
	  if (match(hit)) {
	    top[hit.getModuleID()].insert(hit.getPMTAddress());
	  }
	}
	
	struct JHit_t {
	  
	  double getZ() const { return z; }
	  double getP() const { return p; }
	  
	  double z;
	  double p;
	};
	
	vector<JHit_t> data;
	
	for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
	  
	  if (summary.hasSummaryFrame(module->getID())) {

	    const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
	    
	    JPosition3D pos(module->getPosition());
	    
	    pos.transform(R, tz.getPosition());
	    
	    if (pos.getX() <= roadWidth_m) {
	      
	      const double z = pos.getZ()  -  pos.getX() / getTanThetaC();
	      
	      const double p = getProbability(*module, frame, rates_Hz.getMultiplesRates(), JRegressor_t::T_ns.getLength(), top[module->getID()]);
	      
	      data.push_back({ z, p });
	    }
	  }
	}
	
	double Zmin       = 0.0;     // relative start position [m]
	double Zmax       = 0.0;     // relative end   position [m]
	double npe_total  = 0.0;     // npe due to mip
	double npe_missed = 0.0;     // npe due to mip not detected
	
	if (!data.empty()) {
	  
	  sort(data.begin(), data.end(), JLANG::make_comparator(&JHit_t::getZ));
	  
	  vector<JHit_t>::const_iterator         track_start = start.find(data. begin(), data. end());
	  vector<JHit_t>::const_reverse_iterator track_end   = start.find(data.rbegin(), data.rend());
	  
	  if (track_start != data. end()) { Zmin = track_start->getZ(); }
	  if (track_end   != data.rend()) { Zmax = track_end  ->getZ(); }
	  
	  // additional features
	  
	  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
	    
	    JPosition3D pos(module->getPosition());
	    
	    pos.transform(R, tz.getPosition());
	    
	    if (pos.getX() <= roadWidth_m) {
	      
	      const double z = pos.getZ()  -  pos.getX() / getTanThetaC();
	      
	      if (z >= Zmin && z <= Zmax) {
		
		for (size_t i = 0; i != module->size(); ++i) {
		  
		  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
		  
		  if (getDAQStatus(frame, *module, i)      &&
		      getPMTStatus(frame, *module, i)      &&
		      !module->getPMT(i).has(PMT_DISABLE)) {
		    
		    JPMT pmt = module->getPMT(i);
		    
		    pmt.transform(R, tz.getPosition());
		    
		    const double ya = this->getNPE(pmt, 0.0).getYA();
		    
		    npe_total += ya;

		    if (top[module->getID()].count(i) == 0) {
		      npe_missed += ya;
		    }
		  }
		}
	      }
	    }
	  }
	}
 
	JFit fit = *track;

	// move track
	  
	fit.move(Zmin, getSpeedOfLight());
	  
	out.push_back(fit.add(JMUONSTART));
	  
	out.rbegin()->setW(track->getW());
	out.rbegin()->setW(JSTART_NPE_MIP_TOTAL,  npe_total);
	out.rbegin()->setW(JSTART_NPE_MIP_MISSED, npe_missed);
	out.rbegin()->setW(JSTART_LENGTH_METRES,  Zmax - Zmin);
      }
      
      return out;
    }

  private:
    const JModuleRouter&       router;
    const JSummaryRouter&      summary;
    const JK40Rates            rates_Hz;
    int                        debug;  
  };
}

#endif