//////////////////////////////////////////////////////////////////////////////
// Class:  RAT::Methods::PMTChgpdf
//
// Brief:  Calculate the pdfs of the recorded PMT charge
//
// Author: W. Heintzelman <billh@hep.upenn.edu>
//
// REVISION HISTORY:
//     29/04/2015 : W. Heintzelman - New file
//
// Detail:
//      Given a PMT's characteristics (high-half point, threshold, and
//      pedestal spread) and an assumed number of photoelectrons (PE) produced
//      in  the PMT, calculate two probabilities -- the probability of observing
//      a particular charge value and the probability that the observed charge
//      value will be less than the threshold.  The method takes advantage of
//      the fact that the characteristic function of a sum of independent
//      random variables is the product of the characteristic functions of the
//      variables, and uses Fast Fourier Transform (FFT) methods to calculate
//      and invert characteristic functions.
//
// Terminology:
//      Generated charge is the charge generated in the PMT.  Recorded charge
//      is the charge recorded, i.e., after application of the threshold and
//      the pedestal spread.
//
// Important note:
//      To avoid "wrap-around pollution", there is a joint constraint on the
//      number of points in the discrete pdf input to the FFT (nPts), the
//      charge step between the points (deltaQ), the largest charge for which a
//      non-zero value of the single-PE pdf can be considered (maxSnglQ), the
//      maximum number of photoelectrons considered (dimPE), and the range of
//      non-zero pedestal spread values.  (See "Treatment of End Effects by
//      Zero Padding" in section 13.1 of Numerical Methods in Fortran.)
//      As currently coded, the routine takes into account non-zero pedestal
//      spread pdf values out to 4*sigma, in effect has maxSnglQ= 200, and
//      allows for up to fMaxPE photoelectrons.  The constraint is that
//
//              nPts > (dimPE*maxSnglQ + 4*maxPedSigma)/deltaQ
//
//////////////////////////////////////////////////////////////////////////////


#ifndef __RAT_PMTChgpdf_
#define __RAT_PMTChgpdf_

#include <TFFTComplexReal.h>
#include <TFFTRealComplex.h>
#include <RAT/PMTCharge.hh>
#include <vector>

const int dimPE=15;  // dimension-1 of storage arrays (max allowable no. of PEs)
const int nPts=4096;    // no. of input points to discrete Fourier transforms
const int nOv2P1 = nPts/2 +1;

// Storage for discrete points of a pdf
struct DiscretePDF{
    double p[nPts];
};

// Storage for Fourier coefficients of a discrete characteristic function
struct CharFn{
    double real[nOv2P1];
    double imag[nOv2P1];
};

namespace RAT {

class PMTChgpdf{
    public:
        static PMTChgpdf* Instance();   // Return pointer to singleton
        double PdfRecordedVal
            ( int nPE, double q, double hhp, double threshold, double pedSig) ;
                            // Probability of observing the charge q if nPE
                            // photoelectrons were produced
        double ProbBelowThreshold ( int nPE, double hhp, double threshold) ;
                        // Probability the observed charge is below threshold
        double DoublePolyaQSpectrum(const double Q, double hhp) ;

        void CalcPDFs( double hhp, double threshold, double pedSig) ;
        void CalcPDFs( double hhp, double threshold ) ;

        double DPQSIntegrand (double* x, double* p) ;
        int SetMaxPE(int _maxPE);
        int GetDimPE() { return dimPE; }
        int GetDistnNPts() { return nPts;}
        double GetChgDelta() { return deltaQ;}
        // Warning:  Before calling GetGenChgPDF, the function
        // CalcPDFs(hhp, threshold)  must have been called with the desired
        // values of hhp and threshold.
        double* GetGenChgPDF( int nPE )  { return genChgPDF[nPE]->p ;}

    private:
        static PMTChgpdf* ptrSelf; // Ptr to singleton of this class
        PMTChgpdf() ;              // private because this is a singleton class
        int CheckNPE(int nPE);     // checks that number of PE requested is
                                   // within range provided for
        TFFTRealComplex* tfftFwd;  // FFT from pdf to characteristic fn.
        TFFTComplexReal* tfftRev;  // FFT from char. fn. to pdf
        CharFn* genChgCharFn[dimPE+1];  // characteristic fns of generated chg
        CharFn  spreadCharFn;           // char fn of pedestal spread
        CharFn* truncChgCharFn[dimPE+1];  // char fns of truncated charge,
                                          // before applying pedestal spread
        // In the following pdfs, the value xxxChgPDF[nPE]->p[i] is the
        // probability of observing a charge in the interval of width deltaQ
        // centered on the charge value i*deltaQ, given nPE produced
        // photoelectrons.
        DiscretePDF* genChgPDF[dimPE+1];   // pdfs of generated charge
        DiscretePDF* recChgPDF[dimPE+1];   // pdfs of recorded charge
        PMTCharge* fPMTChargePtr;  // pointer to class having a function to
                        // calculate the single-PE generated-charge probability

        void SetSinglePENorm(double hhp); // Normalize the single-PE charge pdf
        void MultiplyCharFn( CharFn& a, CharFn& b, CharFn& result);
        double PdfBinSum( double* y, int lowLim, int hiLim) ;
        int fMaxPE;      //  Max number of PEs for which calculations are done
        double fnorm;   // normalization of the single PE distn  for last
                        // high half-point
        double hhpLast;  // Last value of high half-point used
        double hhpLastNorm;  // Last value of high half-point used in
                             // normalizing the single-PE charge pdf
        double thresholdLast;  // Last value of threshold
        double spreadLast;  // Last value of pedestal spread
        double deltaQ;      // charge increment for steps in discrete pdfs
};

}   // namespace RAT

#endif