#ifndef __JDETECTOR__JCLBDEFAULTSIMULATOR__
#define __JDETECTOR__JCLBDEFAULTSIMULATOR__

#include <cmath>

#include "JDetector/JCLBDefaultSimulatorInterface.hh"

#include "JTools/JPolint.hh"
#include "JTools/JElement.hh"
#include "JTools/JGridCollection.hh"
#include "JTools/JCollection.hh"
#include "JTools/JGrid.hh"


/**
 * \author mdejong
 */

namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }

namespace JDETECTOR {

  /**
   * Auxiliary class for a non-linear transfer function of TDC inside FPGA.
   */
  template<class JAbscissa_t, class JOrdinate_t>
  class JTransferFunction1D :
    public JTOOLS::JPolintFunction1D<0, 
				     JTOOLS::JElement2D<JAbscissa_t, JOrdinate_t>,
				     JTOOLS::JGridCollection,
				     JOrdinate_t>
  {
  public:

    typedef JTOOLS::JPolintFunction1D<0, 
				      JTOOLS::JElement2D<JAbscissa_t, JOrdinate_t>,
				      JTOOLS::JGridCollection,
				      JOrdinate_t>  JFunction1D_t;


    /**
     * Default constructor.
     */
    JTransferFunction1D()
    {}


    /**
     * Set the non-linearity function.
     *
     * \param  nx              number of elements
     * \param  Xmin            lower limit
     * \param  Xmax            upper limit
     * \param  __begin         begin of weights
     * \param  __end           end   of weights
     */
    template<class T>
    void set(const int         nx,
	     const JAbscissa_t Xmin,
	     const JAbscissa_t Xmax,
	     T __begin,
	     T __end)
    {
      // make CDF

      const int N = std::distance(__begin, __end);

      std::vector<JAbscissa_t> X;

      for (T i = __begin; i != __end; ++i) {
	X.push_back(*i);
      }
      
      for (int i = 0, j = 1; j != N; ++i, ++j) {
	X[j] += X[i];
      }
      
      for (int i = 0; i != N; ++i) {
	X[i] /= X[N-1];
      }


      JFunction1D_t::configure(JTOOLS::make_grid(nx, Xmin, Xmax));

      for (typename JFunction1D_t::iterator i = this->begin(); i != this->end(); ++i) {

	int j = 0;

	while (j != N && Xmin + X[j] * (Xmax - Xmin) < i->getX()) { 
	  ++j;
	}

	i->getY() = (JOrdinate_t) (Xmin  +  j * (Xmax - Xmin) / N);
      }
   
      JFunction1D_t::compile();
    }


    /**
     * Get the function value including non-linearity.
     *
     * \param  x               abscissa value
     * \return                 function value
     */
    JOrdinate_t operator()(const JAbscissa_t x) const
    {
      const long long int n  = (long long int) (x / (this->getXmax() - this->getXmin()));
      const JAbscissa_t   x1 = x  -  n * (this->getXmax() - this->getXmin());

      return (JOrdinate_t) (x - x1)  +  JFunction1D_t::operator()(x1);
    }
  };


  /**
   * Default CLB simulation.
   *
   * This class provides for the implementation of the conversion of a time-over-threshold pulse
   * to a so-called hit.\n
   * A hit consists of:
   *   - PMT readout channel;
   *   - time stamp of the leading edge of the time-over-threshold pulse;
   *   - length of the time-over-threshold pulse;
   * 
   * The time stamping and pulse length determination is done by a TDC 
   * which is integrated inside the FPGA.
   * The observed non-linearity of the TDC is implemented in this class.
   *
   * Although the hits from the same PMT are striclty time sorted,
   * hits from different PMTs could be reordered due to the occupancy of the buffers inside the FPGA.
   * This reordering is also implemented in the this class.
   */
  class JCLBDefaultSimulator :
    public JCLBDefaultSimulatorInterface
  {
  public:
    /**
     * Implementation of non-linearity of TDC.
     */
    class JTDC :
      public JCLBDefaultSimulatorInterface::JTDC
    {
    public:
      /**
       * Default constructor.
       * The non-linearity parameters are due to D. Calvo.
       */
      JTDC()
      {
	const double Tmin_ns = 0.0;
	const double Tmax_ns = 4.0;

	const double W[] = { 0.2413,
			     0.2492,
			     0.2554,
			     0.2541 };

	getTDC.set(4000, Tmin_ns, Tmax_ns, W, W + sizeof(W)/sizeof(W[0]));
      }
      

      /**
       * Make DAQ hit.
       *
       * \param  pmt            PMT channel
       * \param  t_ns           time of hit         [ns]
       * \param  tot_ns         time over threshold [ns]
       * \return                DAQ hit       
       */
      virtual JDAQHit makeHit(const JPMT_t pmt, 
			      const double t_ns,
			      const JTOT_t tot_ns) const override
      {
	return JDAQHit(pmt, getTDC(t_ns), tot_ns);
      }


      JTransferFunction1D<double, JTDC_t> getTDC;     //!< TDC non-linearity function
    };


    /**
     * Auxiliary class to mimic hit ordering effects due to state machine inside CLB.
     */
    class JStateMachine :
      public JCLBDefaultSimulatorInterface::JStateMachine
    {
    public:
      /**
       * Constructor.
       *
       * \param Tmax         maximal time difference for swap [ns]
       */
      JStateMachine(const JDAQHit::JTDC_t Tmax)
      {
	this->Tmax = Tmax;
      }

  
      /**
       * Test whether two consecutive hits may be swapped.
       *
       * \param  first       first  DAQ hit
       * \param  second      second DAQ hit
       * \return             true if PMTs differ and time difference less than preset value; else false
       */
      virtual bool maybeSwapped(const JDAQHit& first, const JDAQHit& second) const override 
      {
	return (first.getPMT() != second.getPMT() && second.getT() - first.getT() < Tmax);
      }


      JDAQHit::JTDC_t Tmax;
    };


    /**
     * Constructor.
     *
     * \param Tmax         maximal time difference for swap [ns]
     */
    JCLBDefaultSimulator(const JDAQHit::JTDC_t Tmax = 0) : 
      JCLBDefaultSimulatorInterface(new JTDC(), 
				    new JStateMachine(Tmax))
    {}
  };
}

#endif