#include <RAT/HitTimeProbabilityFactor.hh>

namespace RAT
{
namespace Methods
{

HitTimeProbabilityFactor::HitTimeProbabilityFactor( int maxPE ) {
    fMaxPE = maxPE;
    fQpdf.resize(fMaxPE+1);
    fFactorial.resize(fMaxPE+1);
    fFactorial[0] = 1.;
    for(int i=1; i<=fMaxPE; i++) {
        fFactorial[i] = (double)i*fFactorial[i-1];
    }
}

void HitTimeProbabilityFactor::SetChargeParameters (
                    double threshold, double highHalfPt, double charge) {
    fPMTChgPdfPtr = RAT::PMTChgpdf::Instance();
    fPMTChgPdfPtr->CalcPDFs(highHalfPt, threshold);
    for(int i=1; i<=fMaxPE; i++) {
        fQpdf[i] =fPMTChgPdfPtr->GetGenChgPDF(i);
                                            // generated charge pdf for i PEs
    }
    double deltaQ = fPMTChgPdfPtr->GetChgDelta();
    fThrIdx = (threshold/deltaQ)+0.5;         // threshold index
    fQIdx = (charge/deltaQ)+0.5;         //  charge index
}

void HitTimeProbabilityFactor::SetTimeParameters(
                                     double pdfPtVal, double cdfPtVal) {
    fPdfPtVal = pdfPtVal; // Value of the time pdf at the observed hit time
    fCdfPtVal = cdfPtVal; // Value of the time cdf at the observed hit time
}

double HitTimeProbabilityFactor::GetTimeProbabilityFactor( int nPE ) {
    // For a given number of PEs, loop over the PEs in the order of their
    // arrival, calculating for each one the probability that it caused the PMT
    // discriminator to fire (i.e., that with the charge it produced  the
    // accumulated charge exceeded the PMT threshold) given the observed
    // charge, times the probability of the observed hit time given that the
    // discriminator fired upon its arrival.  Accumulate the sum of the
    // probabilty products, which is then the probability of the observed hit
    // time given the oberved charge and assuming that fNPE photoelectrons were
    // generated.
    // Functions SetChargeParameters and SetTimeParameters must
    // have been called before calling this function.

    fNPE = nPE;
    double psum = 0.;
    double pFactor= fPdfPtVal*fFactorial[fNPE]*pow(1.-fCdfPtVal,fNPE)/fCdfPtVal;
    double f = fCdfPtVal/(1.-fCdfPtVal);
    for(int iPE=1; iPE<=fNPE; iPE++) {
        // Probability that the iPE-th PMT caused the discriminator to fire
        double pFire = GetPEFireProb( iPE);
        pFactor = f*pFactor;
        psum += pFire*pFactor/(fFactorial[iPE-1]*fFactorial[fNPE-iPE]);
        // The value of the above expression beginning with "pFactor" is the
        // probability of the observed time given that the iPE-th PE caused the
        // discriminator to fire.
    }
    return psum;
}

double HitTimeProbabilityFactor::GetPEFireProb(int firePE) {
    // Given PMT charge (fQIdx), PMT threshold (fThrIdx), the pdf for a PE's
    // generated charge (pdf) already calculated given the PMT's high-half
    // point, and a number of PEs (fNPE) assumed to have been generated,
    // calculate the probabilty that the discriminator fired on the fFirePE-th
    // photoelectron out of fNPE photoelectrons.
    // See SNO+-doc-2871 for an explanation of the computation.
    fFirePE = firePE;
    double prob=0.;
    if(fNPE==1) {
        prob = 1.;
    }else{
        if(fFirePE==1) {
            prob = PSumF(0);
        }else if(fFirePE<fNPE) {    // qp = sum of chgs of PEs preceding fFirePE
            for(int qp=1; qp<=fThrIdx-1; qp++) {
                prob += PQSum(fFirePE-1,qp)*PSumF(qp);
            }
        }else{  //  fFirePE = fNPE
            for(int qp=1; qp<=fThrIdx-1; qp++) {
                prob += PQSum(fFirePE-1,qp)*PQSum(1,fQIdx-qp);
            }
        }
        if( PQSum(fNPE,fQIdx)>0.) prob = prob/PQSum(fNPE,fQIdx);
    }
    return prob;
}

double HitTimeProbabilityFactor::PSumF(int qp) {
    double psum = 0.;       // qf - chg of PE that caused pmt to fire (fFirePE)
    for(int qf=fThrIdx-qp; qf<=fQIdx-qp-(fNPE-fFirePE); qf++) {
        psum += PQSum(1,qf)*PQSum(fNPE-fFirePE, fQIdx-qp-qf);
    }
    return psum;
}

double HitTimeProbabilityFactor::PQSum(int PECount, int q) {
    // return the probability of total charge q from PECount PEs.
    double* probArray = fQpdf[PECount];
    return probArray[q];
}


}       // namespace Methods
}       // namespace RAT