#ifndef __JDETECTOR__JPMTPARAMETERSTOOLKIT__
#define __JDETECTOR__JPMTPARAMETERSTOOLKIT__

#include <limits>

#include "JDetector/JPMTParameters.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JLang/JException.hh"
#include "JLang/JThrow.hh"


/**
 * \author mdejong
 * \file
 * Auxiliary methods to convert hit probability to QE and vice versa.
 */

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

namespace JDETECTOR {

  using JLANG::JValueOutOfRange;


  /**
   * Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assuming QE = 1.
   *
   * \param  parameters        PMT parameters
   * \return                   probability
   */
  inline double getSurvivalProbability(const JPMTParameters& parameters)
  {
    JPMTParameters data = parameters;
    const int      NPE  = 1;

    data.QE = 1.0;

    return JPMTAnalogueSignalProcessor(data).getSurvivalProbability(NPE);
  }


  /**
   * Get ratio of hit probabilities for given QE and expectation value of the number of photo-electrons.\n
   * The ratio corresponds to the hit probability with given QE divided by that with QE = 1.\n
   * The expectation value is defined for a two-fold (or higher) coincidence.
   *
   * \param  QE                QE
   * \param  mu                expectation value
   * \return                   ratio
   */
  inline double getHitProbability(const double QE, const double mu)
  {
    if (mu > 0.0)
      return (1.0 - exp(-QE*mu)) / (1.0 - exp(-mu));
    else
      return QE;
  }


  /**
   * Get maximal ratio of hit probabilities for given QE and expectation value of the number of photo-electrons.\n
   * The ratio corresponds to the hit probability with infinite QE divided by that with QE = 1.\n
   * The expectation value is defined for a two-fold (or higher) coincidence.
   *
   * \param  mu                expectation value
   * \return                   ratio
   */
  inline double getMaximalHitProbability(const double mu)
  {
    using namespace std;

    if (mu > 0.0)
      return 1.0 / (1.0 - exp(-mu));
    else
      return numeric_limits<double>::max();
  }


  /**
   * Get QE for given ratio of hit probabilities and expectation value of the number of photo-electrons.\n
   * The ratio corresponds to the hit probability with given QE divided by that with QE = 1.\n
   * The expectation value is defined for a two-fold (or higher) coincidence.
   *
   * \param  R                 ratio
   * \param  mu                expectation value
   * \return                   QE
   */
  inline double getQE(const double R, const double mu)
  {
    if (R < getMaximalHitProbability(mu)) {
      
      if (mu > 0.0)
	return -log(1.0 - R*(1.0 - exp(-mu))) / mu;
      else
	return R;

    } else {

      THROW(JValueOutOfRange, "Inconsistent arguments at method getQE(" << R << "," << mu <<")");
    }
  }
}

#endif