#include <RAT/DB.hh>
#include <RAT/Log.hh>
#include <RAT/PMTTransitTime.hh>
#include <RAT/PMTTimeDistns.hh>

using std::vector;

namespace RAT {

PMTTimeDistns* PMTTimeDistns::fPtrSelf=NULL;

PMTTimeDistns::PMTTimeDistns() {}

void PMTTimeDistns::BeginOfRun() {
    fPrevPMTIdx = 1000000;
    DB* db = DB::Get();
    DBLinkPtr fLdaq = db->GetLink( "DAQ" );
    fTriggerWindow = fLdaq->GetD( "triggergate" );
    fTriggerToFECDelay = fLdaq->GetD("gtriggerdelay");

    fftFwd = new TFFTRealComplex(fNPts,false);
    fftFwd->Init("M", 0, 0);

    fftRev = new TFFTComplexReal(fNPts,false);
    fftRev->Init("M", 0, 0);

    PMTTransitTime* pmtTransitTimePtr = PMTTransitTime::Get();

    DiscretePDFT pdfTransTime;  // pmtTransTime[i] will be the probability of
        // a PMT transit time in the interval (i-0.5)*fBinWidth to
        // (i+0.5)*fBinWidth.  For a negative time with corresponding
        // (negative) value i, the probability is stored in
        // pmtTransTime[fNPts+i].
    DiscretePDFT pdfBucketTime;  // pmtBucketTime will be the probability of
        // the time from entering the concentrator to creation of a PE being
        // in the interval (i-0.5)*fBinWidth to (i+0.5)*fBinWidth
    CharFnT transTimeCharFn;  // characteristic fn of PMT transit time
    CharFnT bucketTimeCharFn; // characteristic fn of bucket time

    // Calculate discrete pdf of the PMT transit time distribution
    double tSmallest, tLargest;
    pmtTransitTimePtr->GetTimeLimits( tSmallest, tLargest);
    int iMin = floor( tSmallest/fBinWidth );
    int iMax = ceil ( tLargest/fBinWidth );

    for(int i=0; i<fNPts; i++) { pdfTransTime.p[i] =  0.; }
    double hwidth = 0.5*fBinWidth;
    for(int i=iMin; i<iMax; i++) {
        double midpoint = (float)i*fBinWidth;
        int iStore = i;               // Storage for negative times is at
        if (i<0) iStore = fNPts+i;    // upper end of array
        pdfTransTime.p[iStore] =  pmtTransitTimePtr->GetIntervalProbability
                                           ( midpoint-hwidth, midpoint+hwidth );
    }

    // Calculate the discrete characteristic fn of the PMT transit-time pdf:
    fftFwd->SetPoints(pdfTransTime.p);
    fftFwd->Transform();
    fftFwd->GetPointsComplex(transTimeCharFn.real, transTimeCharFn.imag,false);

    // pdf of bucket transit times at 0.2 ns intervals, first bin centered on
    // zero, out to bin centered at 7.8 ns.   This could be put into
    // the database and loaded from there, but there's no advantage to doing so
    // as it does not depend on the run and is not expected to change.
    const int nbPts = 40;
    double bucketpdf[nbPts] = {
        0, 0, 0.395585, 0.468426, 0.0224048, 0.0303111, 0.0188315, 0.0287093,
        0.0107609, 0.00562686, 0.00505185, 0.00535989, 0.00262861, 0.00170449,
        0.000575008, 0.000492864, 0.00115002, 0.00100626, 0.000431256,
        0.000123216, 0.000143752, 0.00020536, 0.000123216, 8.2144e-05,
        2.0536e-05, 4.1072e-05, 6.1608e-05, 0, 0, 2.0536e-05, 0, 0,
        2.0536e-05, 2.0536e-05, 0, 0, 0, 6.1608e-05, 2.0536e-05, 0};

    for(int i=0; i<nbPts; i++) { pdfBucketTime.p[i] = bucketpdf[i]; }
    for(int i=nbPts; i<fNPts; i++) { pdfBucketTime.p[i] = 0.; }

    // Calculate the discrete characteristic fn of the bucket transit-time pdf:
    fftFwd->SetPoints(pdfBucketTime.p);
    fftFwd->Transform();
    fftFwd->GetPointsComplex(bucketTimeCharFn.real,
                                            bucketTimeCharFn.imag,false);

    // Multiply to get characteristic fn. of sum of PMT transit time and bucket
    // transit time
    MultiplyCharFn( transTimeCharFn, bucketTimeCharFn, fDelayCharFn);

}       // PMTTimeDistns::PMTTimeDistns

// *************************************************************************

double PMTTimeDistns::GetRelHitTimeProb( double nExpNoisePEs,
            vector<double>& tCross, vector<double>& tWeight,
            double relHitTime, UInt_t pmtIdx ) {
    // relHitTime is the time of the hit relative to the fitted event time
    if(pmtIdx != fPrevPMTIdx) {
        CalculateTimeDFs( nExpNoisePEs, tCross, tWeight);
        fPrevPMTIdx = pmtIdx;
    }
    // pdf bin in which the relative hit time in the event being fitted falls:
    double adjTime = relHitTime + fTimeOffset;
    int iBin = adjTime/fBinWidth;
    if(iBin<0) iBin = 0;
    if(iBin>(int)fLastBin-1) iBin = (int)fLastBin-1;
    return fHitTimePDF.p[iBin];
}
// *************************************************************************

double PMTTimeDistns::GetRelHitTimeCDF( double nExpNoisePEs,
            vector<double>& tCross, vector<double>& tWeight,
            double relHitTime, UInt_t pmtIdx ) {
    if(pmtIdx != fPrevPMTIdx) {
        CalculateTimeDFs( nExpNoisePEs, tCross, tWeight);
        fPrevPMTIdx = pmtIdx;
    }
    double cdfPtVal;
    double adjTime = relHitTime + fTimeOffset;
    int iBin = adjTime/fBinWidth;
    if(iBin<0) iBin = 0;
    if(iBin>(int)fLastBin-1) {
        iBin = (int)fLastBin-1;
        cdfPtVal = fHitTimeCDF.p[iBin];
    }else{
        // Calculate the value of the cumulative distribution function at the
        // lower and upper edge of ibin (dividing pdf[] by fNPts is necessary
        // to get proper normalization of the FFT output pdf):
        cdfPtVal = fHitTimeCDF.p[iBin-1] +
                   fHitTimePDF.p[iBin]*(adjTime - iBin*fBinWidth)/fBinWidth;
    }
    return cdfPtVal;
}
// *************************************************************************


void PMTTimeDistns::CalcOrderStats( double nExpNoisePEs,
    vector<double>& tCross, vector<double>& tWeight,
    double relHitTime, UInt_t pmtIdx,
    double pWidth, double chgIntT, double tLimit, vector< vector<double> >& g,
    vector< vector<double> >& Gwid, vector< vector<double> >& Gprm,
    vector< vector<double> >& Glim) {
    // Calculate the time distributions of the order statistics of
    // the PE times.
    // pWidth is the pulse width
    // chgIntT is the length of the charge integration time interval.
    // tLimit is the estimated latest time for a hit (not a PE) to be included
    //          in the event
    // fMaxPE is the maximum number of PEs for which the calculation is done.
    // Results:  g[j][i] = probability that the time of the ith out of j PEs
    //                     is in the time bin corresponding to relHitTime.
    //           Gwid[j][i] = probability that the time of the ith out of j PEs
    //                     is less than or equal to (relHitTime + pWidth),
    //                     given that the first one is at relHitTime.
    //           Gprm[j][i] = probability that the time of the ith out of j PEs
    //                     is less than or equal to (relHitTime + chgIntT),
    //                     given that the first one is at relHitTime.
    //           Glim[j][i] = probability that the time of the ith out of j PEs
    //                     is less than or equal to tLimit
    if(pmtIdx != fPrevPMTIdx) {
        CalculateTimeDFs( nExpNoisePEs, tCross, tWeight);
        fPrevPMTIdx = pmtIdx;
    }

    int widthBin = (relHitTime + pWidth + fTimeOffset)/fBinWidth;
    int lastChgBin  = (relHitTime + chgIntT + fTimeOffset)/fBinWidth;
    int hitBin   = (relHitTime          + fTimeOffset)/fBinWidth;
    int limitBin = (tLimit              + fTimeOffset)/fBinWidth;
    int lastBin = widthBin;
    if(lastChgBin>lastBin) lastBin = lastChgBin;
    if(hitBin>lastBin) lastBin = hitBin;
    if(limitBin>lastBin) lastBin = limitBin;

    vector< vector<double> > FP, CP;
    // FP[i][n] becomes the integral, F[i],  of the single-PE time pdf f[i] up
    // to the time corresponding to bin i, raised to the nth power.
    // CP[i][n] is (1-F[i]) raised to the nth power.
    FP.resize( lastBin+1, vector<double>(fMaxPE+1,0.) );
    CP.resize( lastBin+1, vector<double>(fMaxPE+1,0.) );

    // tidx is the time index
    for(int tidx=0; tidx<=lastBin; tidx++) {
        double F = fHitTimeCDF.p[tidx];
        double C = 1.-F;
        if(F==0.) {
            FP[tidx][0] = 0.;
            CP[tidx][0] = 1.;
        }else if (F==1.){
            FP[tidx][0] = 1.;
            CP[tidx][0] = 0.;
        }else{
            FP[tidx][0] = 1.;
            CP[tidx][0] = 1.;
        }
        double Fprod = 1.;
        double Cprod = 1.;
        for(int kpow=1; kpow<=fMaxPE; kpow++) {
            Fprod = Fprod*F;
            Cprod = Cprod*C;
            FP[tidx][kpow] = Fprod;
            CP[tidx][kpow] = Cprod;
        }
    }

    Gwid.clear();
    Gwid.resize(fMaxPE+1);
    Gprm.clear();
    Gprm.resize(fMaxPE+1);
    Glim.clear();
    Glim.resize(fMaxPE+1);
    g.clear();
    g.resize(fMaxPE+1);
    for(int j=1; j<=fMaxPE; j++) {
        g[j].push_back(0.);        // g[j][0] = 0.;
        Gwid[j].push_back(0.);     // Gwid[j][0] = 0.;
        Gprm[j].push_back(0.);     // Gprm[j][0] = 0.;
        Glim[j].push_back(0.);     // Glim[j][0] = 0.;

        for(int i=1; i<=j; i++) {
            // Integrate order statistic pdfs over the time bins:
            double sum=0.;
            double GTHit = 0.;
            double GTWidth, GTChg;
            for(int tidx=0; tidx<=lastBin; tidx++) {
                double term = FP[tidx][i-1]*CP[tidx][j-i]*fHitTimePDF.p[tidx];
                sum = sum + term;
                if (tidx==hitBin) {
                    double value = i*fBinomCoeff[j][i]*term ;
                    g[j].push_back( value );
                    GTHit = i*fBinomCoeff[j][i]*sum ;
                }
                if (tidx==widthBin)  GTWidth = i*fBinomCoeff[j][i]*sum ;
                if (tidx==lastChgBin) GTChg = i*fBinomCoeff[j][i]*sum;
                if (tidx==limitBin) {
                    double value =  i*fBinomCoeff[j][i]*sum;
                    if(value>1.) value = 1.;
                    Glim[j].push_back( value);
                }
            }

            double Gwidth;
            if(GTHit >= GTWidth ) {
                Gwidth = 0.;
            }else if(GTHit >= 1. ) {
                Gwidth = 1.;
            }else{
                Gwidth = (GTWidth-GTHit)/(1.-GTHit);
            }

            double Gprime;
            if(GTHit >= GTChg ) {
                Gprime = 0.;
            }else if(GTHit >= 1. ) {
                Gprime = 1.;
            }else{
                Gprime = (GTChg-GTHit)/(1.-GTHit);
            }

            if(Gwidth<0.) info << "Gwidth<0 " <<" "<<Gwidth <<" "<< j <<" "<< i
                <<" "<< GTHit <<" "<< sum <<" "<< GTChg << endl;
            if(Gwidth>1.) Gwidth = 1.;    // cumulative errors can cause sum>1
            Gwid[j].push_back( Gwidth);         // Gwid[j][i] =

            if(Gprime<0.) info << "Gprime<0 " <<" "<<Gprime <<" "<< j <<" "<< i
                <<" "<< GTHit <<" "<< sum <<" "<< GTChg << endl;
            if(Gprime>1.) Gprime = 1.;
            Gprm[j].push_back( Gprime);         // Gprm[j][i] =
            if(pmtIdx==2717) {
            }

            if(Glim[j][i]>1.) Glim[j][i] = 1.;
        }       // for(int i=1; i<=j; i++)
        Gwid[j][1] = 1.;
        Gprm[j][1] = 1.;
        Gwid[j].push_back( 0.); // Gwid[j][j+1] = 0. for calling routine to work
        Gprm[j].push_back( 0.); // Gprm[j][j+1] = 0. for calling routine to work
    }
    return;
}

// *************************************************************************

void PMTTimeDistns::CalcIntervalProbs(double nExpNoisePEs,
                    vector<double>& tCross, vector<double>& tWeight,
                    double relHitTime, UInt_t pmtIdx,
                    double delta, vector< vector<double> >& pdelta) {
    // pdelta[n][i] is the probability, given that n PEs occur after
    // relHitTime, that i of those will occur before relHitTime+delta
    if(pmtIdx != fPrevPMTIdx) {
        CalculateTimeDFs( nExpNoisePEs, tCross, tWeight);
        fPrevPMTIdx = pmtIdx;
    }

    double FtH = GetRelHitTimeCDF( nExpNoisePEs, tCross, tWeight,
                                                    relHitTime, pmtIdx);
    double FtHdelta = GetRelHitTimeCDF( nExpNoisePEs, tCross, tWeight,
                                                    relHitTime+delta, pmtIdx);
    pdelta.clear();
    pdelta.resize(fMaxPE+1);
    pdelta[0].push_back(1.);            // pdelta[0][0] = 1
    if(FtH==1.) {  // All of the time distn from the auxiliary simulation
                   // precedes the recorded hit time.  Shouldn't happen, but we
                   // need to do something reasonable if it does because of low
                   // statistics.
        for(int n=1; n<=fMaxPE; n++) {
            pdelta[n].push_back( 0. );      // pdelta[n][0] = 0
            for(int i=1; i<=n; i++) {
                pdelta[n].push_back( 1. );  // pdelta[n][i] = 1, for i>=1
            }
        }
    }else{
        double eta = (FtHdelta - FtH)/(1. - FtH);
        for(int n=1; n<=fMaxPE; n++) {
            double factor = 1.;
            for(int i=0; i<=n; i++) {
                pdelta[n].push_back(fBinomCoeff[n][i]*factor); // pdelta[n][i] =
                factor *= eta;
            }
            factor = 1.;
            for(int i=n; i>=0; i--) {
                pdelta[n][i] *= factor;
                factor *= (1. - eta);
            }
        }
    }
    return;
}

// *************************************************************************

void PMTTimeDistns::CalcBinomCoeff(int maxn, vector< vector<double> >& c) {
// Calculate array of binomial coefficients up to maximum value of n
    c.clear();
    c.resize(maxn+1);

    c[0].push_back(1.);                          //  c[0][0] = 1;
    for(int n=1; n<=maxn; n++) {
        c[n].push_back(1.);                      //  c[n][0] = 1;
        for(int i=1; i<=n; i++) {
            // c[n][i] = c[n][i-1]*(double)(n-i+1)/(double)i;
            c[n].push_back( c[n][i-1]*(n-i+1)/i );
        }
    }
    return;
}

// *************************************************************************

void PMTTimeDistns::CalculateTimeDFs( double nExpNoisePEs,
        vector<double>& tCross, vector<double>& tWeight) {
    // Note:  The value of fEstTriggerTime for the event must have been set by
    // a call to SetEventParameters before this function is called.

    // Calculate the pdf of the crossing times from the array tCross, which is
    // a list of the crossing times.  Use tWeight to weight the times.
    // Add possible noise hits to the pdf.
    DiscretePDFT pdfCrossTime;
    for(int i=0; i<fNPts; i++) { pdfCrossTime.p[i] = 0.; }

    size_t nCross = tCross.size();
    if (nCross != tWeight.size()) warn << "PMTTimeDistns::CalculateTimeDFs: "
        << "Size inconsistency -- "<< nCross <<" "<< tWeight.size() << endl;
    double totalWeight = 0.;
    for(size_t i=0; i<nCross; i++) {
        totalWeight += tWeight[i];
    }

    // Hack -- scale nExpNoisePEs to length of noise window used here:
    // 400 is triggergate -- should get this number out of the database instead
    // of hardcoding it here.
    nExpNoisePEs = nExpNoisePEs*fNoiseLimit/400.;

    // Fraction of PEs created in PMT expected to be from noise (which is then
    // equal to the contribution of noise PEs to the pdf of the hit time
    // distribution):
    double fNoise = nExpNoisePEs/(nExpNoisePEs + totalWeight );
    int firstNoiseBin = fTimeOffset/fBinWidth;
    int lastNoiseBin = (fNoiseLimit + fTimeOffset)/fBinWidth;
    if (lastNoiseBin>3700) {
        info << "lastNoiseBin larger than expected" << lastNoiseBin<< endl;
        lastNoiseBin=3700;
    }

    // Distribute total noise probability evenly across bins in the trigger
    // window.  The noise hit probability is scaled to sum to the correct total
    // within the bins corresponding to the trigger window.
    double baseProb = fNoise/ (lastNoiseBin-firstNoiseBin);
    for(int i=firstNoiseBin; i<lastNoiseBin; i++) {
        pdfCrossTime.p[i] = baseProb;
    }

    if(nCross != 0) {
        // Distribute total probability fraction for signal evenly among bins
        // where an incident photon was observed in the auxiliary simulation.
        // We depend on the transit-time spread to smooth the probability
        // over multiple bins.
        double f = (1.-fNoise)/totalWeight;
        for(size_t  i=0; i<nCross; i++) {
            int ibin = (tCross[i] + fTimeOffset)/fBinWidth;
            pdfCrossTime.p[ibin] += f*tWeight[i];
        }
    }
    // Calculate the discrete characteristic fn of the crossing time pdf:
    CharFnT crossTimeCharFn;
    fftFwd->SetPoints(pdfCrossTime.p);
    fftFwd->Transform();
    fftFwd->GetPointsComplex(crossTimeCharFn.real,
                                            crossTimeCharFn.imag,false);
    // Calculate the discrete characteristic fn of the total time pdf:
    CharFnT totalTimeCharFn;
    MultiplyCharFn( fDelayCharFn, crossTimeCharFn, totalTimeCharFn);

    // Inverse FFT to get pdf of total time:
    fftRev->SetPointsComplex( totalTimeCharFn.real, totalTimeCharFn.imag );
    fftRev->Transform();
    double* pdf = fftRev->GetPointsReal(false);
    double fInvNPts = 1./(double)fNPts;

    // In the following, avoid possible small negative pdf values that
    // could be produced by the FFT.
    // Dividing by fNPts is necessary to get proper normalization of
    // the FFT output pdf.
    if( pdf[0]>0. ) {
        fHitTimePDF.p[0] = pdf[0]*fInvNPts;
        fHitTimeCDF.p[0] = fHitTimePDF.p[0];
    }else{
        fHitTimePDF.p[0] = 1.e-33;
        fHitTimeCDF.p[0] = 0.;
    }
    for(int j=1; j<fNPts; j++) {
        if( pdf[j]>0. ) {
            fHitTimePDF.p[j] = pdf[j]*fInvNPts;
            fHitTimeCDF.p[j] = fHitTimeCDF.p[j-1] +fHitTimePDF.p[j];
        }else{
            fHitTimePDF.p[j] = 1.e-33;
            fHitTimeCDF.p[j] = fHitTimeCDF.p[j-1];
        }
        // Ensure that cumulative errors don't cause CDF to exceed 1.
        if( fHitTimeCDF.p[j] > 1.) fHitTimeCDF.p[j] = 1.;
    }
}

// *************************************************************************

void PMTTimeDistns::SetEventParameters ( double estTriggerTime, size_t maxPE ){
    // estTriggerTime is the estimated event trigger time relative to event
    // start time.
    fEstTriggerTime = estTriggerTime;
//     fTimeOffset = fTriggerWindow - (fEstTriggerTime + fTriggerToFECDelay );
    fHighLimit = fEstTriggerTime + fTriggerToFECDelay + 390. -10.;
    // The long integration time is currently hardcoded as 390 ns, but 10 ns of
    // that is before the PMT fires.
    fNoiseLimit = fHighLimit;
    fLastBin = fHighLimit/fBinWidth;
    fTimeOffset = 40.;
    fMaxPE = maxPE;
    CalcBinomCoeff( fMaxPE, fBinomCoeff );
}

// *************************************************************************

PMTTimeDistns* PMTTimeDistns::Instance() {
    if(fPtrSelf==NULL) { fPtrSelf = new PMTTimeDistns(); }
    return fPtrSelf;
}

// *************************************************************************

void PMTTimeDistns::MultiplyCharFn( CharFnT& a, CharFnT& b, CharFnT& result) {
    // Multiplication of two characteristic functions
    for(int i=0; i<fNOv2P1; i++) {
        result.real[i] = a.real[i]*b.real[i] - a.imag[i]*b.imag[i];
        result.imag[i] = a.imag[i]*b.real[i] + a.real[i]*b.imag[i];
    }
}

}   // namespace RAT