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

using namespace std;

namespace RAT {

PMTTimeDistns* PMTTimeDistns::fPtrSelf=NULL;

PMTTimeDistns::PMTTimeDistns() {}

void PMTTimeDistns::BeginOfRun() {
    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

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

void PMTTimeDistns::GetRelHitTimeProbs(  double nExpNoisePEs, double relHitTime,
            vector<double>& tCross, double& pdfPtVal, double& cdfPtVal ) {
    CalculateValues( nExpNoisePEs, relHitTime,  tCross) ;
    pdfPtVal = fPdfPtVal;
    cdfPtVal = fCdfPtVal;
}
// *************************************************************************

void PMTTimeDistns::GetRelHitTimeCDF(double nExpNoisePEs, double relHitTime,
            vector<double>& tCross, double& cdfA, double& cdfB ) {
    CalculateValues( nExpNoisePEs, relHitTime,  tCross) ;
    cdfA = fCdfA;
    cdfB = fCdfB;
}

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

void PMTTimeDistns::GetRelHitTimeCCDF(double nExpNoisePEs, double relHitTime,
            vector<double>& tCross, double& ccdfA, double& ccdfB ) {
    CalculateValues( nExpNoisePEs, relHitTime,  tCross) ;
    ccdfA = 1. - fCdfA;
    ccdfB = 1. - fCdfB;
}

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

void PMTTimeDistns::CalculateValues( double nExpNoisePEs, double relHitTime,
            vector<double>& tCross) {
    // Note:  The values of fEstTriggerTime, fLowCut, and fHighCut 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.  Add possible noise hits to the pdf.
    DiscretePDFT pdfCrossTime;
    for(int i=0; i<fNPts; i++) { pdfCrossTime.p[i] = 0.; }

    int nWindowBins = ceil((fTriggerWindow)/fBinWidth);  // bins spanned by
                                                         // trigger window
    size_t nCross = tCross.size();

    // 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 + (double)nCross );
    int firstBin = (fLowCut + fTimeOffset)/fBinWidth;
    if (firstBin<0) firstBin=0;
    int lastBin = (fHighCut + fTimeOffset)/fBinWidth;
    if (lastBin>nWindowBins) lastBin=nWindowBins;
    // Distribute total noise probability evenly across bins in the trigger
    // window, or if the ModeCut selector is used, across the reduced interval.
    // The noise hit probability is scaled to sum to the correct total within
    // the bins corresponding to the trigger window or reduced interval.
    double baseProb = fNoise/ (lastBin-firstBin);
    for(int i=firstBin; i<lastBin; i++) { pdfCrossTime.p[i] = baseProb; }

//     detail << "PMTTimeDistns: fEstTriggerTime, fTimeOffset, nWindowBins, " <<
//            "fNoise, firstBin, lastBin = " <<fEstTriggerTime<<" "<<
//                fTimeOffset<<" "<<  nWindowBins<<" "<<  fNoise<<" "<<
//                firstBin<<" "<<  lastBin << endl;

    if(nCross != 0) {
        // Distribute total probability fraction for signal evenly among bins
        // where an incident photon was observed in the auxiliary simulation.
        double f = (1.-fNoise)/(double)nCross;
        for(size_t  i=0; i<nCross; i++) {
            int ibin = (tCross[i] + fTimeOffset)/fBinWidth;
            pdfCrossTime.p[ibin] += f;
        }
    }
    // 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;

    // 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>=nWindowBins-1) iBin = nWindowBins-2;
    // 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):
    double sumProbLower = 0.;
    for(int j=0; j<iBin; j++) {
        // In the following, avoid possible small negative pdf values that
        // could be produced by the FFT.
        if( pdf[j]>0. ) sumProbLower += pdf[j];
    }
    sumProbLower *= fInvNPts;
    double sumProbUpper = sumProbLower;
    fPdfPtVal =  0.;
    if( pdf[iBin]>0.) {
        fPdfPtVal =  pdf[iBin]*fInvNPts;
        sumProbUpper += fPdfPtVal;
    }
    fCdfA = sumProbLower;
    fCdfB = sumProbUpper;
    fCdfPtVal = sumProbLower + fPdfPtVal*(adjTime - iBin*fBinWidth)/fBinWidth;
}       // void PMTTimeDistns::CalculateValues

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

void PMTTimeDistns::SetEventParameters
                    ( double estTriggerTime, double lowCut, double highCut ) {
    // estTriggerTime is the estimated trigger time relative to event start
    // time.
    fEstTriggerTime = estTriggerTime;
    fLowCut = lowCut;
    fHighCut = highCut;
    fTimeOffset = fTriggerWindow - (fEstTriggerTime + fTriggerToFECDelay );
//     detail <<"SetEventParameters: estTriggertime, fTriggerToFECDelay, " <<
//         "fTimeOffset,lowCut, highCut = " <<
//          estTriggerTime<<" "<< fTriggerToFECDelay<<" "<< fTimeOffset<<" "<<
//          lowCut<<" "<<  highCut << endl;
}

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

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