//////////////////////////////////////////////////////////////////////////////
// Class:  RAT::Methods::PMTTimeDistns
//
// Brief:  Calculate various probabilities related to an observed PMT hit time,
// for 0 to fMaxPE possible generated photoelectrons
//
// Author: W. Heintzelman <billh@hep.upenn.edu>
//
// REVISION HISTORY:
//     25/05/2016 : W. Heintzelman - New file
//
// Detail:
//      Given an input table that defines the pdf of the time a photon spends
//      in the PMT bucket (i.e., after crossing the PMT face) until a PE is
//      produced, the distribution of the PMT transit time (i.e., the time
//      interval between when a PE is produced and the recorded hit time) from
//      class PMTTransitTime, and a list of tabulation-area photon-crossing
//      times from the auxiliary simulation, this class computes values of the
//      complementary cumulative distribution of the total time between the
//      occurence of an event and the recorded hit time.  It uses 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.
//
// Notes: The times corresponding to values in the time pdfs are relative to
//      the estimated beginning of the trigger window of the event, which is
//      calculated from an estimate of the trigger time.
//
//      This is a singleton class
//
//////////////////////////////////////////////////////////////////////////////


#ifndef __RAT_PMTTimeDistns_
#define __RAT_PMTTimeDistns_

#include <TFFTComplexReal.h>
#include <TFFTRealComplex.h>
#include <vector>

using std::vector;

const int fNPts=4096;    // no. of input points to discrete Fourier transforms
const int fNOv2P1 = fNPts/2 +1;
const double fBinWidth = 0.2;  // bin width in ns of input and output time pdfs

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

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

namespace RAT {

class PMTTimeDistns{
    public:
        static PMTTimeDistns* Instance();   // Return pointer to singleton
        void   BeginOfRun();

        void SetEventParameters
                        (double estTriggerTime, double lowCut, double highCut);

        // Following function has as inputs:
        //    nExpNoisePEs  - The auxiliary simulation is executed without
        //           noise.  This is the expected number of noise PEs in the
        //           PMT in question during the trigger window if noise had
        //           been modelled in the auxiliary simulation, multiplied by
        //           the expected ratio of the number of photons incident on
        //           PMT tabulation area to the number of photoelectrons
        //           produced by photons.  This ratio is calculated in PEnergy
        //           as equal to the tabulation area factor times the number of
        //           initiating tracks in the auxiliary simulation divided by
        //           the average PMT effectiveness factor for the incident
        //           photons.
        //    relHitTime = the hit time of the PMT relative to fitted event time
        //    tCross = list of times photons crossed the tabulation area
        //              of the PMT in the auxiliary simulation, relative to
        //              the start of the MC event.
        //
        // Outputs:

        //     cdfPtVal  = value of the cumulative distribution function at the
        //                      hit time
        //     pdfPtVal  = value of the probability density
        //           function at the hit time
        void GetRelHitTimeProbs(  double nExpNoisePEs, double relHitTime,
                vector<double>& tCross, double& pdfPtVal, double& cdfPtVal );

        //     cdfA  = value of the cumulative distribution
        //           function of the hit time, at the beginning of the time bin
        //           containing the recorded PMT hit time.
        //     cdfB  = value of the CDF at the end of the time bin.
        void GetRelHitTimeCDF(  double nExpNoisePEs, double relHitTime,
                vector<double>& tCross, double& cdfA, double& cdfB );

        //     ccdfA  = value of the complementary cumulative distribution
        //           function of the hit time, at the beginning of the time bin
        //           containing the recorded PMT hit time.
        //     ccdfB  = value of the CCDF at the end of the time bin.
        void GetRelHitTimeCCDF(  double nExpNoisePEs, double relHitTime,
                vector<double>& tCross, double& ccdfA, double& ccdfB );

    private:
        static PMTTimeDistns* fPtrSelf; // Ptr to singleton of this class
        PMTTimeDistns() ;          // private because this is a singleton class

        TFFTRealComplex* fftFwd;  // FFT from pdf to characteristic fn.
        TFFTComplexReal* fftRev;  // FFT from char. fn. to pdf
        CharFnT fDelayCharFn;       // characteristic fn of total time between
                // when photon enters concentrator and recorded hit on FEC.

        void CalculateValues( double nExpNoisePEs, double relHitTime,
                                                    vector<double>& tCross) ;
            // The values of fEstTriggerTime, fLowCut, and fHighCut for the
            // event must have been set by a call to SetEventParameters before
            // this function is called.

        void MultiplyCharFn( CharFnT& a, CharFnT& b, CharFnT& result);
        double fTimeOffset;     // Offset applied to photon crossing times
                    // before doing the FFT so time at lower edge of the
                    // trigger window is zero rather than negative.  As a
                    // result, the times corresponding to the bins in the pdfs
                    // are relative  to the lower edge of the trigger window.

        double fTriggerWindow;  // length of the trigger window (gate) in ns
        double fTriggerToFECDelay; // delay of trigger signal back to FEC in ns
        double fEstTriggerTime;  // estimated trigger time of auxiliary event
                                 // relative to the event start time
        double fLowCut, fHighCut;  // Limits of time cut window if ModeCut or
                           // TimeResidualCut is used, otherwise large defaults
        double fCdfA, fCdfB, fPdfPtVal, fCdfPtVal;
};

}   // namespace RAT

#endif