#ifndef __JPHYSICS__JK40RATES___
#define __JPHYSICS__JK40RATES___

#include <vector>

#include "JMath/JMathSupportkit.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {

  /**
   * Type definition of singles rate [Hz].
   */
  typedef double                 JRateL0_t;
  
  /**
   * Type definition of count rate as a function of multiplicty [Hz]
   * The multiples rate start counting at two-fold coincidences.
   */
  typedef std::vector<double>    JRateL1_t;


  /**
   * Type definition of multiplicity.
   */
  typedef size_t                 multiplicity_type;

  
  /**
   * Auxiliary class for K40 rates.
   * The singles rate refers to the counting rate of a PMT and
   * the multiples rate to the rate of genuine coincidences due to K40 decays.
   */
  struct JK40Rates {
    /**
     * Default constructor.
     */
    JK40Rates() :
      rateL0(0.0),
      rateL1()
    {}


    /**
     * Constructor.
     *
     * The multiples rates start counting at two-fold coincidences.
     *
     * \param  rateL0_Hz       singles   rate  [Hz]
     * \param  rateL1_Hz       multiples rates [Hz]
     */
    JK40Rates(const JRateL0_t  rateL0_Hz,
	      const JRateL1_t& rateL1_Hz = JRateL1_t()) :
      rateL0(rateL0_Hz),
      rateL1(rateL1_Hz)
    {}


    /**
     * Get singles rate.
     *
     * \return                 rate [Hz]
     */
    double getSinglesRate() const
    {
      return rateL0;
    }


    /**
     * Get multiples rate.
     *
     * \return                 rate [Hz]
     */
    const JRateL1_t& getMultiplesRates() const
    {
      return rateL1;
    }


    /**
     * Get multiples rate at given multiplicity.
     *
     * \param  M               multiplicity (M >= JK40Rates::LOWER_L1_MULTIPLICITY)
     * \return                 rate [Hz]
     */
    double getMultiplesRate(const multiplicity_type M) const
    {
      if (M >= LOWER_L1_MULTIPLICITY && M - LOWER_L1_MULTIPLICITY < (multiplicity_type) rateL1.size())
	return rateL1[M - LOWER_L1_MULTIPLICITY];
      else
	return 0.0;
    }


    /**
     * Get lower multiplicty.
     *
     * \return                 lower multiplicity
     */
    multiplicity_type getLowerL1Multiplicity() const
    {
      return LOWER_L1_MULTIPLICITY;
    }


    /**
     * Get upper multiplicty.
     *
     * \return                 upper multiplicity
     */
    multiplicity_type getUpperL1Multiplicity() const
    {
      return rateL1.size() + 1;
    }

    
    /**
     * Correct rates for global efficiency,
     *
     * \param  QE              global efficiency
     */
    void correct(const double QE)
    {
      if (QE > 0.0) {
      
	rateL0 /= QE;

	JRateL1_t buffer = rateL1;
	
	for (multiplicity_type M = getUpperL1Multiplicity(); M >= getLowerL1Multiplicity(); --M) {

	  // determine contribution from higher multiplicities
	  
	  double R = 0.0;
	  
	  for (multiplicity_type i = M + 1; i <= getUpperL1Multiplicity(); ++i) {
	    R += buffer[i - LOWER_L1_MULTIPLICITY] * JMATH::binomial(i, M) * pow(QE, M) * pow(1.0 - QE, i - M);
	  }

	  if (getMultiplesRate(M) > R)
	    buffer[M - LOWER_L1_MULTIPLICITY] = (getMultiplesRate(M) - R) / pow(QE, M);
	  else
	    buffer[M - LOWER_L1_MULTIPLICITY] =  0.0;
	}

	rateL1 = buffer;
	
      } else {

	rateL0 = 0.0;

	for (JRateL1_t::iterator i = rateL1.begin(); i != rateL1.end(); ++i) {
	  *i = 0.0;
	}
      }
    }
    
    
    /**
     * Read K40 rates from input.
     *
     * \param  in              input stream
     * \param  object          K40 rates
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JK40Rates& object)
    {
      const double rateL0 = object.rateL0;

      if (in >> object.rateL0) {

	object.rateL1.clear();

	for (double x; in >> x; ) {
	  object.rateL1.push_back(x);
	}

      } else {

	object.rateL0 = rateL0;
      }

      return in;
    }


    /**
     * Write K40 rates to output.
     *
     * \param  out             output stream
     * \param  object          K40 rates
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JK40Rates& object)
    {
      out << object.rateL0;

      for (JRateL1_t::const_iterator i = object.rateL1.begin(); i != object.rateL1.end(); ++i) {
	out << ' ' << *i;
      }

      return out;
    }

    
    /**
     * Lower L1 multiplicity
     */
    static const multiplicity_type  LOWER_L1_MULTIPLICITY  =  2;

    
  protected:
    JRateL0_t rateL0;          //!< singles   rate  [Hz]
    JRateL1_t rateL1;          //!< multiples rates [Hz]
  };
}

#endif