// ChanSWStatusProc.cc
// Contact person: Billy Liggins <billy.liggins@qmul.ac.uk>
// See ChanSWStatusProc.hh for more details
//———————————————————————//

#include <TFile.h>

#include <RAT/Processor.hh>
#include <RAT/ChanSWStatusProc.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DS/UniversalTime.hh>
#include <RAT/DB.hh>
#include <RAT/DBTable.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/ChanHWStatus.hh>
#include <RAT/DB.hh>
#include <RAT/DBLink.hh>

#include <TVectorD.h>
#include <TTree.h>
#include <TF1.h>

#include <map>
#include <string>
#include <iterator>
#include <ctime>
#include <time.h>
#include <exception>

using namespace std;


namespace RAT
{

    void ChanSWStatusProc::BeginOfRun( DS::Run& run )
    {

        fRunID = run.GetRunID();
        GetStartTime();
        fHighOccEvents = 0;
        chargeBins=5096;
        chargeMin=-1000;
        chargeMax=4096;

        OccBins=1e5;
        OccMin=1e-7;
        OccMax=1.0;
        //30 mins
        fTimeInterval = 30.0 * 60 * 1e9;
        //Returns total number of channels in the database.
        fAllLCNs = DU::Utility::Get()->GetPMTInfo().GetCount();

        // Read in cuts from CSS_CUTS_VALUES.ratdb
        DBLinkPtr lCSSCUTS = DB::Get()->GetLink("CSS_CUT_VALUES", "SNO_default");
        fHighOccCut = lCSSCUTS->GetD( "high_occ" );
        fLowOccCut=  lCSSCUTS->GetD( "low_occ" );

        // Filling the CHS vector with status pass/fail for every channel
        const DU::ChanHWStatus& channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();
        fCHSVector.assign(fAllLCNs, 0); // all online
        for(size_t lcn = 0; lcn < fAllLCNs; lcn++)
        {
            if( !channelHardwareStatus.IsTubeOnline(lcn) ){
                fCHSVector.at(lcn)=1;
                continue;
            }
            // And create a histogram for each channel for QHS, QLX, QHL, TAC
            hQHS[lcn] = new TH1D(dformat("QHS_proc_%u", lcn).c_str(), dformat("QHS_proc_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hQHL[lcn] = new TH1D(dformat("QHL_proc_%u", lcn).c_str(), dformat("QHL_proc_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hQLX[lcn] = new TH1D(dformat("QLX_proc_%u", lcn).c_str(), dformat("QLX_proc_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hTAC[lcn] = new TH1D(dformat("TAC_proc_%u", lcn).c_str(), dformat("TAC_proc_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hGainRatio[lcn] = new TH1D(dformat("GainRatio_proc_%u", lcn).c_str(), dformat("GainRatio_proc_%u", lcn).c_str(), 220,-10,100);
        }
    }

    Processor::Result ChanSWStatusProc::DSEvent( DS::Run&, DS::Entry& ds )
    {
        for( size_t iEV = 0; iEV < ds.GetEVCount(); iEV++ ){
            const bool pedEvent = BitManip::TestBit( static_cast<unsigned int>( ds.GetEV(iEV).GetTrigType() ),
                    DU::Utility::Get()->GetTrigBits().GetBitIndex( "Pedestal" ) );
            if(!pedEvent)
                Event(ds, ds.GetEV(iEV));
        }
        return OK;
    }

    Processor::Result ChanSWStatusProc::Event( DS::Entry&, DS::EV& ev )
    {
        const bool pulseGTEvent = BitManip::TestBit( static_cast<unsigned int>( ev.GetTrigType() ),
                DU::Utility::Get()->GetTrigBits().GetBitIndex( "PulseGT" ) );

        // Counting the total number of PGT events
        if( pulseGTEvent )
            fHighOccEvents++;

        const DS::UniversalTime delta = ev.GetUniversalTime() - fStartRun;
        const int interval = delta.GetNanoSeconds() / fTimeInterval;

        // If fLowOccEvents doesn't have index interval,
        // make one and make sure that all lower intervals are initialised.
        if ( (size_t) interval >= fLowOccEvents.size() ){
            fLowOccEvents.resize(interval+1);
        }

        // Counts the number of events in each interval.
        fLowOccEvents.at(interval)++;

        // The ECA and PCA status are stored in the same bit mask. The first 32 bits
        // are dedicated to the statuses of ECA checks, while the following 32 bits
        // contain the statuses of PCA checks
        RAT::DS::BitMask PCAStatus = DS::BitMask();
        RAT::DS::BitMask ECAStatus = DS::BitMask();

        for( size_t ipmt = 0; ipmt < ev.GetCalPMTs().GetCount(); ipmt++ )
        {
            const DS::PMTCal& pmt = ev.GetCalPMTs().GetPMT(ipmt);
            unsigned int lcn = pmt.GetLCN();

            // Offline channels are excluded from the study
            if( fCHSVector.at(lcn) == 1 ) continue;

            // Filling the empty HighOccMap
            // This is where Jeanne wants to change the way te low occ is calulated.
            if( pulseGTEvent )
            {
                if( fHighOccMap.find( lcn ) == fHighOccMap.end() )
                    fHighOccMap[lcn] = 0; // If lcn was not found, then create an entry
                fHighOccMap[lcn]++; // Number of hits from PGT events for each PMT
            }

            // Filling the empty LowOccMap if involved in an event.
            if( fLowOccMap[interval].find( lcn ) == fLowOccMap[interval].end() )
                fLowOccMap[interval][lcn] = 0; // If this lcn not yet initialised for this interval, do so
            fLowOccMap[interval][lcn]++; // Number of hits from events for each PMT

            // Filling histograms with charge and time for each event per LCN
            hQHS[lcn]->Fill( pmt.GetQHS() );
            hQHL[lcn]->Fill( pmt.GetQHL() );
            hQLX[lcn]->Fill( pmt.GetQLX() );
            hTAC[lcn]->Fill( pmt.GetTime() );
            hGainRatio[lcn]->Fill( pmt.GetQHL() / pmt.GetQLX() );

            // Filling empty PCA map. indexed by lcn and containing another map, that
            // contains strings of statuses and number of hits correspoding to them
            PCAStatus = pmt.GetStatus();
            if( fPCAMap[lcn].find( "total" ) == fPCAMap[lcn].end() )
            {
                fPCAMap[lcn]["total"] = 0;
                fPCAMap[lcn]["failed"] = 0;
            }
            fPCAMap[lcn]["total"]++;

            // Bit 32 is resposible for overall PCA status
            // If is is equal to 1, that means it failed
            if( PCAStatus.Get(32) == 1)
                fPCAMap[lcn]["failed"]++;
        }

        //This loop is to fill the ECA map.

        for( size_t ipmt = 0; ipmt < ev.GetCalPMTs().GetAllCount(); ipmt++ )
        {
            const DS::PMTCal& pmt = ev.GetCalPMTs().GetAllPMT(ipmt);
            unsigned int lcn = pmt.GetLCN();
            ECAStatus = pmt.GetStatus();

            // Offline channels are excluded from the study
            if( fCHSVector.at(lcn) == 1 ) continue;

            // Filling empty ECA map. indexed by lcn and containing another map, that
            // contains strings of statuses and number of hits correspoding to them
            if( fECAMap[lcn].find( "total" ) == fECAMap[lcn].end() )
            {
                fECAMap[lcn]["total"] = 0;
                fECAMap[lcn]["failed"] = 0;
            }
            fECAMap[lcn]["total"]++;

            // Bit 0 is resposible for overall ECA status
            // If is is equal to 1, that means it failed
            if( ECAStatus.Get(0) == 1)
                fECAMap[lcn]["failed"]++;
        }
        return OK;
    }

    void ChanSWStatusProc::GetStartTime(){
        RAT::DB *db = RAT::DB::Get();
        RAT::DBLinkPtr dblink = db->GetLink("RUN");

        const int day = dblink->GetI("start_day");
        const int secs= dblink->GetI("start_sec");
        const int nano = static_cast<int>( dblink->GetD("start_nsc"));
        fStartRun = DS::UniversalTime(day,secs,nano);
    }

    void ChanSWStatusProc::CalculateECACut( TFile& outputFile, RAT::DBTable& table ) const
    {
        std::vector<int> failed(fAllLCNs, DS::INVALID);
        std::vector<double> ecaFailedFraction(fAllLCNs, DS::INVALID);

        for( CalibLCNCountMap::const_iterator iLCN = fECAMap.begin(); iLCN != fECAMap.end(); iLCN++ )
        {
            failed.at(iLCN->first) = (iLCN->second).at("failed");
            ecaFailedFraction.at(iLCN->first) = static_cast<double>( (iLCN->second).at("failed") ) /static_cast<double>( (iLCN->second).at("total") );
        }
        table.SetDArray("eca_failed_fraction", ecaFailedFraction);
        table.SetIArray("eca_failed", failed);
    }


    void ChanSWStatusProc::CalculatePCACut( TFile& outputFile, RAT::DBTable& table ) const
    {
        std::vector<int> failed(fAllLCNs, DS::INVALID);
        std::vector<double> pcaFailedFraction(fAllLCNs, DS::INVALID);

        for( CalibLCNCountMap::const_iterator iLCN = fPCAMap.begin(); iLCN != fPCAMap.end(); iLCN++ )
        {
            failed.at(iLCN->first) =  (iLCN->second).at("failed");
            pcaFailedFraction.at(iLCN->first) = static_cast<double>( (iLCN->second).at("failed") ) /static_cast<double>( (iLCN->second).at("total") );
        }
        table.SetDArray("pca_failed_fraction", pcaFailedFraction);
        table.SetIArray("pca_failed", failed);
    }

    double ChanSWStatusProc::GetFWHM(const TH1D& hist_) const{
      const int bin1 = hist_.FindFirstBinAbove(hist_.GetMaximum()/2);
      const int bin2 = hist_.FindLastBinAbove(hist_.GetMaximum()/2);
      return hist_.GetBinCenter(bin2) - hist_.GetBinCenter(bin1);
    }

    void ChanSWStatusProc::CalculateHighOccCut( TFile& outputFile, RAT::DBTable& table ) const
    {
        TH1D highOccRate("HighOccRate","Channel PGT hits/Number of Pulsed Global triggered events",OccBins,OccMin,OccMax);
        for( LCNCountMap::const_iterator iTer = fHighOccMap.begin(); iTer != fHighOccMap.end(); iTer++ )
        {
            // rate = (number of hits of PGT events of a channel) / (total number of hits of PGT events)
            const double rate = static_cast<double>((*iTer).second) / static_cast<double>(fHighOccEvents);
            highOccRate.Fill(rate);
        }
        outputFile.cd();
        highOccRate.Write();
        const double highOcc_FWHM = GetFWHM(highOccRate);
        const double cutValue = highOccRate.GetMean() + (fHighOccCut * highOcc_FWHM);

        table.SetD( "high_occ_cut", cutValue );
    }


    void ChanSWStatusProc::CalculateLowOccCut( TFile& outputFile, RAT::DBTable& table ) const
    {
        TH1D lowOccRate("LowOccRate","Channel hits per time interval/Number of total hits per time interval",OccBins,OccMin,OccMax);

        for( IntervalLCNCountMap::const_iterator iInterval = fLowOccMap.begin(); iInterval != fLowOccMap.end(); iInterval++ )
        {
            for( LCNCountMap::const_iterator iLCN = iInterval->second.begin(); iLCN != iInterval->second.end(); iLCN++ )
            {
                // rate = (number of hits of a channel per time block) / (total number of events per interval)
                const double rate = static_cast<double>( iLCN->second ) / static_cast<double>(fLowOccEvents.at(iInterval->first));
                lowOccRate.Fill(rate);
            }
        }
        outputFile.cd();
        lowOccRate.Write();

        // Finding the mode of the histogram.
        const Double_t mode = lowOccRate.GetBinCenter(lowOccRate.GetMaximumBin());
        lowOccRate.Fit("gaus","S","",0,mode);
        const TF1 *fit = lowOccRate.GetFunction("gaus");
        const Double_t rms = fit->GetParameter(2);
        const Double_t mean = fit->GetParameter(1);
        const double cutValue = mean - (fLowOccCut * rms);
        table.SetD( "low_occ_cut", cutValue );
    }


    void ChanSWStatusProc::CalculateChargeTimeCut( TFile& outputFile, RAT::DBTable& table )
    {

        // In order to output the cut values into ratdb table in a string temporary vectors
        // filled with DS::INVALID have been created. The cut values of each check are presented
        // by two numbers. Strings of low and high numbers of each check are outputed into the ratdb
        // table separately. If the channel is off, its cut value stays [-9999, -9999], while
        // the cut values of working channels are updated. Since LCN are counted from 1 and vector
        // elements are counted from 0, the first element of every vector is either -9999 or 9999
        // as a default.
        std::vector<double> qhlLowVector(fAllLCNs, DS::INVALID);
        std::vector<double> qhlHighVector(fAllLCNs, DS::INVALID);
        std::vector<double> qhsLowVector(fAllLCNs, DS::INVALID);
        std::vector<double> qhsHighVector(fAllLCNs, DS::INVALID);
        std::vector<double> qlxLowVector(fAllLCNs, DS::INVALID);
        std::vector<double> qlxHighVector(fAllLCNs, DS::INVALID);
        std::vector<double> tacLowVector(fAllLCNs, DS::INVALID);
        std::vector<double> tacHighVector(fAllLCNs, DS::INVALID);

        // Calculate cut values for charge and time and save histograms for each channel
        outputFile.cd();
        for (std::map<unsigned int, TH1D*>::const_iterator ilcn = hQHS.begin(); ilcn !=hQHS.end(); ++ilcn) {

            hQHS[ilcn->first]->Write();
            hQHL[ilcn->first]->Write();
            hQLX[ilcn->first]->Write();
            hTAC[ilcn->first]->Write();
            hGainRatio[ilcn->first]->Write();


            double quantiles[] = { 0.05, 0.95 };
            double cutQHS[] = { DS::INVALID, DS::INVALID };
            double cutQHL[] = { DS::INVALID, DS::INVALID };
            double cutQLX[] = { DS::INVALID, DS::INVALID };
            double cutTAC[] = { DS::INVALID, DS::INVALID };

            // The function GetQuantiles() calculates the low
            // and high cuts using the cumulative sum
            // Returns low and high cuts stored as an array
            if(!( hQHS[ilcn->first]->Integral()==0 or hQHS[ilcn->first]->ComputeIntegral()==0 or
                        hQHS[ilcn->first]->GetEntries()==0))
                hQHS[ilcn->first]->GetQuantiles( 2, cutQHS, quantiles );

            if(!( hQHL[ilcn->first]->Integral()==0 or hQHL[ilcn->first]->ComputeIntegral()==0 or
                        hQHL[ilcn->first]->GetEntries()==0))
                hQHL[ilcn->first]->GetQuantiles( 2, cutQHL, quantiles );

            if(!( hQLX[ilcn->first]->Integral()==0 or hQLX[ilcn->first]->ComputeIntegral()==0 or
                        hQLX[ilcn->first]->GetEntries()==0))
                hQLX[ilcn->first]->GetQuantiles( 2, cutQLX, quantiles );

            if(!( hTAC[ilcn->first]->Integral()==0 or hTAC[ilcn->first]->ComputeIntegral()==0 or
                        hTAC[ilcn->first]->GetEntries()==0))
                hTAC[ilcn->first]->GetQuantiles( 2, cutTAC, quantiles );

            qhlLowVector. at(ilcn->first) = cutQHL[0];
            qhlHighVector.at(ilcn->first) = cutQHL[1];
            qhsLowVector. at(ilcn->first) = cutQHS[0];
            qhsHighVector.at(ilcn->first) = cutQHS[1];
            qlxLowVector. at(ilcn->first) = cutQLX[0];
            qlxHighVector.at(ilcn->first) = cutQLX[1];
            tacLowVector. at(ilcn->first) = cutTAC[0];
            tacHighVector.at(ilcn->first) = cutTAC[1];
        }
        table.SetDArray( "qhl_low", qhlLowVector );
        table.SetDArray( "qhl_high", qhlHighVector );
        table.SetDArray( "qhs_low", qhsLowVector );
        table.SetDArray( "qhs_high", qhsHighVector );
        table.SetDArray( "qlx_low", qlxLowVector );
        table.SetDArray( "qlx_high", qlxHighVector );
        table.SetDArray( "tac_low", tacLowVector );
        table.SetDArray( "tac_high", tacHighVector );
    }


    void ChanSWStatusProc::EndOfRun( DS::Run& run )
    {
        TFile outputFile( dformat("CSSStandardRun_%d.root",run.GetRunID()).c_str(), "recreate");
        RAT::DBTable table( "CSS_CUTS" );
        table.SetI( "version", 2 );
        table.SetRunRange( fRunID, RAT::DB::DBINFINITY ); // online mode table. Valid until infinity
        table.SetPassNumber( -2 ); // online mode table. Pass is assigned by the database server
        table.SetS( "comment", dformat( "Created by ChanSWStatusProc applied to %d", fRunID ) );
        table.SetIArray("chs_status", fCHSVector);
        CalculateHighOccCut( outputFile, table );
        CalculateLowOccCut( outputFile, table );
        CalculateChargeTimeCut( outputFile, table );
        CalculateECACut( outputFile, table );
        CalculatePCACut( outputFile, table );
        // From https://stackoverflow.com/questions/9527960/how-do-i-construct-an-iso-8601-datetime-in-c
        time_t now;
        time( &now );
        char timestamp[sizeof "2015-07-02T14:27:06Z"];
        strftime( timestamp, sizeof timestamp, "%FT%TZ", gmtime( &now ) );
        table.SetS( "timestamp", timestamp );
        table.SaveAs( dformat("CSS_CUTS_%u.ratdb", run.GetRunID()) );
        outputFile.Close();

        // Now tidy up the histograms we created
        for (std::map<unsigned int, TH1D*>::const_iterator lcn = hQHS.begin(); lcn != hQHS.end(); ++lcn) {
            delete hQHS.at(lcn->first); // delete the histogram
            delete hQHL.at(lcn->first);
            delete hQLX.at(lcn->first);
            delete hTAC.at(lcn->first);
            delete hGainRatio.at(lcn->first);
        }

        hQHS.clear();
        hQHL.clear();
        hQLX.clear();
        hTAC.clear();
        hGainRatio.clear();

    }

}