#ifndef __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
#define __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__

#include <vector>
#include <set>

#include "TRandom3.h"

#include "JDetector/JK40Simulator.hh"
#include "JDetector/JModuleIdentifier.hh"
#include "JDetector/JPMTIdentifier.hh"
#include "JDetector/JTimeRange.hh"

#include "JMath/JMathToolkit.hh"
#include "JLang/JComparator.hh"
#include "JLang/JException.hh"
#include "JLang/JCC.hh"


/**
 * \author mdejong
 */

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

namespace JDETECTOR {

  using JLANG::JNumericalPrecision;


  /**
   * Default K40 simulator interface.
   *
   * This class provides for a default implementation of the JK40Simulator interface
   * which is based on a set of virtual methods.
   * These methods constitute a user interface to the K40 background simulation.
   */
  class JK40DefaultSimulatorInterface :
    public JK40Simulator
  {
    /**
     * PMT pair.
     */
    struct pair_type {
      size_t pmt1;     //!< first  PMT
      size_t pmt2;     //!< second PMT
      double P;        //!< probability of coincidence
    };


    /**
     * Auxiliary data structure for the probabilities of genuine coincidences due to K40 decays as a function of the PMT pair.\n
     * The first element in this list should correspond to an arbitray PMT pair with zero probability before using the function operator.
     */
    mutable struct :
      public std::vector<pair_type>
    {
      /**
       * Get random PMT pair according probabilities of genuine coincidences due to K40 decays.
       *
       * \param  random         random value <0,1]
       * \return                PMT pair
       */
      const pair_type& operator()(const double random) const
      {
	double P = this->rbegin()->P;

	if (P <= 0.0) {
	  THROW(JNumericalPrecision, "Probability " << P);
	}

	P *= random;

	int imin = 0;
	int imax = this->size() - 1;

	for (int i = (imax + imin) / 2; imax - imin != 1; i = (imax + imin) / 2) {
	  if ((*this)[i].P < P)
	    imin = i;
	  else
	    imax = i;
	}

	return (*this)[imax];
      }
    } probabilityL1;
    

  protected:
    /**
     * Default constructor.
     */
    JK40DefaultSimulatorInterface()
    {}


  public:
    /**
     * Get singles rate as a function of PMT.
     *
     * \param  pmt            PMT identifier
     * \return                rate [Hz]
     */
    virtual double getSinglesRate(const JPMTIdentifier& pmt) const = 0;

    
    /**
     * Get multiples rate as a function of optical module.
     *
     * \param  module         optical module identifier
     * \param  M              multiplicity (M >= 2)
     * \return                rate [Hz]
     */
    virtual double getMultiplesRate(const JModuleIdentifier& module, const int M) const = 0;


    /**
     * Get probability of coincidence.
     *
     * \param  ct             cosine space angle between PMT axes
     * \return                probability
     */
    virtual double getProbability(const double ct) const = 0;


    /**
     * Generate hits.
     *
     * \param  module         module
     * \param  period         time window [ns]
     * \param  output         background data
     */
    virtual void generateHits(const JModule&    module,
			      const JTimeRange& period,
			      JModuleData&      output) const
    {
      using namespace std;
      using namespace JPP;
      

      // resize internal buffers
      
      const size_t N = module.size();
      const size_t M = (N * (N - 1)) / 2;

      rateL1_Hz    .resize(N);
      probabilityL1.resize(M + 1);


      // generate singles

      for (size_t pmt = 0; pmt != N; ++pmt) {

	const double rateL0_Hz = getSinglesRate(JPMTIdentifier(module.getID(), pmt));

	if (rateL0_Hz > 0.0) {
	      
	  const double t_ns = 1.0e9 / rateL0_Hz;          // [ns]
	      
	  for (double t1 = period.getLowerLimit() + gRandom->Exp(t_ns); t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
	    output[pmt].push_back(JPMTSignal(t1, 1));
	  }
	}
      }


      // generate coincidences

      double totalRateL1_Hz = 0.0;

      for (size_t i = 0; i != N; ++i) {
	totalRateL1_Hz += rateL1_Hz[i] = getMultiplesRate(module.getID(), i + 2);
      }

      if (totalRateL1_Hz > 0.0) {

	const double t_ns = 1.0e9 / totalRateL1_Hz;       // [ns]

	double t1 = period.getLowerLimit() + gRandom->Exp(t_ns);

	if (t1 < period.getUpperLimit()) {

	  // configure pair-wise propabilities

	  probabilityL1[0] = { N, N, 0.0};

	  size_t i = 0;
	  double P = 0.0;

	  for (size_t pmt1 = 0; pmt1 != N; ++pmt1) {
	    for (size_t pmt2 = 0; pmt2 != pmt1; ++pmt2) {

	      const double ct = getDot(module[pmt1].getDirection(), module[pmt2].getDirection());
	      const double p  = getProbability(ct);

	      i += 1;
	      P += p;

	      probabilityL1[i] = { pmt1, pmt2, P};
	    }
	  }
	  
	  for ( ; t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {

	    try {

	      // generate two-fold coincidence

	      const pair_type& pair = probabilityL1(gRandom->Rndm());

	      output[pair.pmt1].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
	      output[pair.pmt2].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
	      
	      // generate larger than two-fold coincidences, if any

	      size_t M = 0;

	      for (double R = totalRateL1_Hz * gRandom->Rndm(); M != N && (R -= rateL1_Hz[M]) > 0.0; ++M) {}

	      if (M != 0) {

		set<size_t> buffer = { pair.pmt1, pair.pmt2 };   // hit PMTs

		for ( ; M != 0; --M) {

		  vector<double> probability1D(N, 1.0);

		  double P = 0.0;
		  
		  for (size_t i = 0; i != N; ++i) {
		    
		    if (buffer.count(i) == 0) {

		      for (set<size_t>::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {

			const double ct = getDot(module[i].getDirection(), module[*pmt].getDirection());

			probability1D[i] *= getProbability(ct);
		      }

		      P += probability1D[i];

		    } else {

		      probability1D[i] = 0.0;
		    }
		  }

		  if (P > 0.0) {

		    size_t pmt = 0;

		    for (P *= gRandom->Rndm(); pmt != N && (P -= probability1D[pmt]) > 0.0; ++pmt) {}

		    if (pmt != N) {

		      output[pmt].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));

		      buffer.insert(pmt);
		    }

		  } else {

		    break;
		  }
		}
	      }
	    }
	    catch (const JNumericalPrecision&) {}
	  }
	}
      }
    }


    /**
     * Get intrinsic time smearing of K40 coincidences.
     *
     * \return                sigma [ns]
     */
    static double getSigma()
    {
      return get_sigma();
    }


    /**
     * Set intrinsic time smearing of K40 coincidences.
     *
     * \param  sigma          sigma [ns]
     */
    static void setSigma(const double sigma)
    {
      get_sigma() = sigma;
    }

  private:
    /**
     * Get intrinsic time smearing of K40 coincidences.
     *
     * \return                sigma [ns]
     */
    static double& get_sigma()
    {
      static double sigma = 0.5;

      return sigma;
    }


    /**
     * Multiples rate as a function of the multiplicity.
     * The index i corresponds to multiplicity M = i + 2.
     */
    mutable std::vector<double> rateL1_Hz;
  }; 
}

#endif