////////////////////////////////////////////////////////////////////
/// \class RAT::PMTCharge
///
/// \brief   Calculate charge induced in a PMT for a single incident photon
///
/// \author Stan Seibert <volsung@physics.utexas.edu>
///
/// REVISION HISTORY:\n
///     ~~ Oct 2008 : Richie Bonventre - Rewrite entire code to include
///                   dipolya charge distribution as done in snoman \n
///     24 Feb 2010 : Gabriel Orebi Gann - Declare PIntegrate public
///                                        Add a few comments to the
///                                        code for clarity
///     10 Sep 2010 : Gabriel Orebi Gann - Impose min/max limits on HHP, check
///     tube status, & set online tubes with bad HHP values to the mean
///     30 Oct 2013 : Matt Mottram - Make a singleton
///     23 Jun 2014 : A. Mastbaum - Add flag to enable/disable channel
///                                 specific HHP, plus check that
///                                 N_hhp == N_pmts
///     13 May 2015 : W Heintzelman - change access of GetHHP to public
///
/// \details The single pe charge distribution is modelled as a
///          `dipolya' plus an exponential. This routine throws a charge
///          from that distribution, to be assigned to a hit PMT. The
///          charge distribution is dependent on the individual tube
///          gains, as determined by the high-half points, read from
///          PMTCALIB.ratdb. Actually throwing a random number (R) from
///          the analytic distribution is slow, so instead a uniform
///          random number is thrown, and the following eqn is solved
///          for Qroot: R - integral[0, Qroot]{P(Q)dQ} = 0.
///
////////////////////////////////////////////////////////////////////
#ifndef __RAT_PMTCharge__
#define __RAT_PMTCharge__

#include <CLHEP/Random/RandGeneral.h>

#include <RAT/DB.hh>

#include <vector>
#include <cstddef>

namespace RAT {

class PMTCharge {
public:

  /// Singleton class instance
  inline static PMTCharge* Get();

  /// Destructor
  virtual ~PMTCharge();

  /// Initialise settings from the ratdb
  void BeginOfRun();

  /// Returns charge for npe photoelectrons
  virtual float CalcCharge(const int npe);

  /// Returns charge for npe photoelectrons fitted to CDF
  virtual float CalcCharge(int npe, bool Polya, const int ipmt);

  /// Returns integral of charge distribution from [0, b] for pmt impt
  double PIntegrate(const double b, const int ipmt);

  /// Average charge spectrum
  float DoublePolyaQSpectrum(const double Qcd );

  /// Per-PMT charge spectrum, dependent on the high half point (read from PMTCALIB.ratdb)
  float DoublePolyaQSpectrum(const double Qcd, double hhp );

  // Get a channel-specific HHP, or the mean if channel-specific values are
  // disabled.
  inline double GetHHP(int i) const;

protected:

  /// \brief Throw a charge according to the PMT single pe charge distribution.
  /// This is done by taking a random number (flat) and finding the root
  /// of the equation which is the difference between flat and the
  /// integral of the charge spectrum up to the root (Method taken from
  /// snoman)
  float Bisection(const double flat, const int ipmt, const double norm);

  /// Return Gamma(A) * (1 - Gamma(X))
  double Gamma(const double A, const double X);

  CLHEP::RandGeneral* fQdist;
  double fQlo, fQhi;
  float fM,  fGamM, fBeta;
  float fA1, fB1, fC1, fQP0, fQPa;
  DBLinkPtr fPmtCalib;
  std::vector<float> fQhshhp;

  // Flag to enable/disable channel-specific HHP values
  bool fUseQHSHHP;

  // The mean HHP value from the database, used if channel-specific is disabled
  double fMeanHHP;
};

inline PMTCharge* PMTCharge::Get()
{
  static PMTCharge instance;
  return &instance;
}

inline double PMTCharge::GetHHP(int i) const
{
  if (fUseQHSHHP) {
    return fQhshhp.at(i);
  }
  return fMeanHHP;
}

} // namespace RAT
#endif