// ChanSWStatusCalib.cc
// Contact person: Lorna Nolan <l.j.nolan@qmul.ac.uk>
// See ChanSWStatusCalib.hh for more details
//———————————————————————//

#include <TH1D.h>
#include <TFile.h>

#include <RAT/Processor.hh>
#include <RAT/ChanSWStatusCalib.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/Log.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/MetaInformation.hh>

#include <map>
#include <string>
#include <iterator>
#include <ctime>

using namespace std;

namespace RAT
{

    void ChanSWStatusCalib::BeginOfRun( DS::Run& run )
    {
        fRunID = run.GetRunID();
        GetStartTime();
        fHighOccEvents = 0;
        chargeBins=5096;
        chargeMin=-1000;
        chargeMax=4096;

        fTimeInterval = 30.0 * 60 * 1e9;
        fInitialised = false;
        fAllLCNs = DU::Utility::Get()->GetPMTInfo().GetCount();
        fCSSBits = DU::Utility::Get()->GetChanSWStatus();

        // Getting the vector of Channel Hardware Status in the standard run
        fCHSVectorStandardRun = DB::Get()->GetLink( "CSS_CUTS" )->GetIArray( "chs_status" );
        // Getting the vector of Channel Hardware Status in the current run
        const DU::ChanHWStatus& channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();
        fCHSVector.assign(fAllLCNs, 0); // all passed
        fChannelWords.assign(fAllLCNs,0); // all checks pass by defualt.
        for(size_t lcn = 0; lcn < fAllLCNs; lcn++)
        {
            if( !channelHardwareStatus.IsTubeOnline(lcn) ){
                fCHSVector.at(lcn)=1; // offline
                continue;
            }
            // And create a histogram for each channel for QHS, QLX, QHL, TAC
            hQHS[lcn] = new TH1D(dformat("QHS_calib_%u", lcn).c_str(), dformat("QHS_calib_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hQHL[lcn] = new TH1D(dformat("QHL_calib_%u", lcn).c_str(), dformat("QHL_calib_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hQLX[lcn] = new TH1D(dformat("QLX_calib_%u", lcn).c_str(), dformat("QLX_calib_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hTAC[lcn] = new TH1D(dformat("TAC_calib_%u", lcn).c_str(), dformat("TAC_calib_%u", lcn).c_str(), chargeBins,chargeMin,chargeMax);
            hGainRatio[lcn] = new TH1D(dformat("GainRatio_calib_%u", lcn).c_str(), dformat("GainRatio_calib_%u", lcn).c_str(), 220,-10,100);
        }
        // get the charge cuts
        DBLinkPtr lCSSCUTS = DB::Get()->GetLink("CSS_CUT_VALUES", "SNO+_default");
        negativeCut = lCSSCUTS->GetD("negative_percent");
        mainCut = lCSSCUTS->GetD("main_percent");
        secondCut = lCSSCUTS->GetD("second_percent");
        thirdCut = lCSSCUTS->GetD("third_percent");
        railedCut = lCSSCUTS->GetD("railed_percent");
        // Fetch the high and low occupancy cuts
        highOccCut = DB::Get()->GetLink( "CSS_CUTS" )->GetD( "high_occ_cut" );
        lowOccCut = DB::Get()->GetLink( "CSS_CUTS" )->GetD( "low_occ_cut" );
        // Fetch the fraction of eca and pca failed
        ecaStandardFF = DB::Get()->GetLink( "CSS_CUTS" )->GetDArray( "eca_failed_fraction" );
        pcaStandardFF = DB::Get()->GetLink( "CSS_CUTS" )->GetDArray( "pca_failed_fraction" );
        // Fetch the charge ranges from table
        negativeChargeRange = lCSSCUTS->GetDArray("qhs_negative_range");
        mainChargeRange = lCSSCUTS->GetDArray("qhs_main_range");
        secondChargeRange = lCSSCUTS->GetDArray("qhs_second_range");
        thirdChargeRange = lCSSCUTS->GetDArray("qhs_third_range");
        railedChargeRange = lCSSCUTS->GetDArray("qhs_railed_range");
    }


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

        // 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();
        // 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 ( static_cast<size_t>(interval) >= fLowOccEvents.size() ){
            fLowOccEvents.resize(interval+1);
        }

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

        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 map
            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( fLowOccMap[interval].find( lcn ) == fLowOccMap[interval].end() )
                fLowOccMap[interval][lcn] = 0;
            fLowOccMap[interval][lcn]++;

            // 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"]++;
        }

        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
    ChanSWStatusCalib::ApplyCHSCut()
        {
            // Setting the chs bit to 1 (default 0) iff the channel fails chs
            // for this run. Setting the chs_standard bit to 1 (default 0) iff
            // the channel failed chs in the standard run.

            for(size_t lcn = 0; lcn < fAllLCNs; lcn++)
            {
                if( fCHSVector.at(lcn)==1 )
                    fChannelWords.at(lcn) |= 1<<fCSSBits.GetBitIndex("chs");

                if( fCHSVectorStandardRun.at(lcn)==1 )
                    fChannelWords.at(lcn) |= 1<<fCSSBits.GetBitIndex("chs_standard");
            }
        }

    void
    ChanSWStatusCalib::ApplyHighOccCut()
        {
            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)
                double rate = static_cast<double>((*iTer).second) / static_cast<double>(fHighOccEvents);
                if( rate > highOccCut )
                    fChannelWords.at(iTer->first) |= 1<<fCSSBits.GetBitIndex("high_occ");
            }
        }

    void
    ChanSWStatusCalib::ApplyECACut()
        {
            std::vector<double> ecaFailedFraction(fAllLCNs, DS::INVALID);

            for( CalibLCNCountMap::const_iterator iLCN = fECAMap.begin(); iLCN != fECAMap.end(); iLCN++ )
            {
                // the ratio between the number of hits that failed ECA over the total number of hits of
                // every channel in a run
                ecaFailedFraction.at(iLCN->first) = static_cast<double>( (iLCN->second).at("failed") ) /static_cast<double>( (iLCN->second).at("total") );

                if( ecaFailedFraction.at(iLCN->first)!= 1.0 and ecaStandardFF.at(iLCN->first)!= 1.0 and
                        abs( ecaFailedFraction.at(iLCN->first) - ecaStandardFF.at(iLCN->first) ) > ( 0.1 * ecaStandardFF.at(iLCN->first) ) ){
                    fChannelWords.at(iLCN->first) |= 1<<fCSSBits.GetBitIndex("eca");
                }
            }
        }

    void
    ChanSWStatusCalib::ApplyPCACut()
        {
            std::vector<double> pcaFailedFraction(fAllLCNs, DS::INVALID);

            for( CalibLCNCountMap::const_iterator iLCN = fPCAMap.begin(); iLCN != fPCAMap.end(); iLCN++ )
            {
                // the ratio between the number of hits that failed PCA over the total number of hits of
                // every channel in a run
                pcaFailedFraction.at(iLCN->first) = static_cast<double>( (iLCN->second).at("failed") ) /static_cast<double>( (iLCN->second).at("total") );
                if( abs( pcaFailedFraction.at(iLCN->first) - pcaStandardFF.at(iLCN->first) ) > ( 0.1 * pcaStandardFF.at(iLCN->first) ) and pcaStandardFF.at(iLCN->first) != 1 and pcaFailedFraction.at(iLCN->first) != 1 ){
                    fChannelWords.at(iLCN->first) |= 1<<fCSSBits.GetBitIndex("pca");
                }

            }
        }


    void
    ChanSWStatusCalib::ApplyLowOccCut()
        {

            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 block)
                    const double rate = static_cast<double>( iLCN->second ) / static_cast<double>(fLowOccEvents.at(iInterval->first));

                    // If a channel fails the check in any time block, the channel is considered to fail in general
                    if( rate < lowOccCut )
                        fChannelWords.at(iLCN->first) |= 1<<fCSSBits.GetBitIndex("low_occ");
                }
            }
        }

    void
    ChanSWStatusCalib::ApplyChargeTimeCut( const  std::map<unsigned int,TH1D*>& histos, const string& type )
        {
            // In order to output the check values into ratdb table, a temporary vector is created.
            // If the channel fails the check, its status is updated.
            std::vector<int> chargeTimeVector(fAllLCNs, 0);
            std::vector<double> cuts_low = DB::Get()->GetLink( "CSS_CUTS" )->GetDArray( dformat("%s_low", type.c_str()) );
            std::vector<double> cuts_high = DB::Get()->GetLink( "CSS_CUTS" )->GetDArray( dformat("%s_high", type.c_str()) );

            //JW Loop through all the channels - one histo for each channel and evaluate fraction that fall outside cuts.
            // If 50% outside cuts for that channel - fail it
            // JW need to check histo integral logic and binning accuracy
            for (std::map<unsigned int, TH1D*>::const_iterator ilcn = histos.begin(); ilcn != histos.end(); ++ilcn) {
                if(cuts_low[ilcn->first] == DS::INVALID || cuts_high[ilcn->first] == DS::INVALID){
                    // What is the correct logic if the cut is invalid? - in previous code ignored and therefore passed
                    // fChannelWords.at(ilcn->first) |= 1<<fCSSBits.GetBitIndex(type);
                    continue;
                }else{
                    double norm = histos.at(ilcn->first)->GetEntries(); // all entries including those outside range in overflows
                    int binlow = histos.at(ilcn->first)->FindBin(cuts_low[ilcn->first]);
                    int binhigh = histos.at(ilcn->first)->FindBin(cuts_high[ilcn->first]);
                    double frac = histos.at(ilcn->first)->Integral(binlow,binhigh)/norm;  // this is the fraction within cuts
                    if(frac<0.5)
                        fChannelWords.at(ilcn->first) |= 1<<fCSSBits.GetBitIndex(type);
                }
            }
        }
    void
    ChanSWStatusCalib::ApplyQHSChargeCut(const  std::map<unsigned int,TH1D*>& histos)
        {
            // In order to output the check values into ratdb table, a temporary vector is created.
            // If the channel fails the check, its status is updated.
            std::vector<int> chargeTimeVector(fAllLCNs, 0);

            //Loop through all channels and then apply cuts
            for (std::map<unsigned int, TH1D*>::const_iterator ilcn = histos.begin(); ilcn != histos.end(); ++ilcn) {
              int b1 = histos.at(ilcn->first)->FindBin(chargeMin);
              int b2 = histos.at(ilcn->first)->FindBin(chargeMax);
              Double_t whole_integral = histos.at(ilcn->first)->Integral(b1, b2);
              int b1Negative = histos.at(ilcn->first)->FindBin(negativeChargeRange[0]);
              int b2Negative = histos.at(ilcn->first)->FindBin(negativeChargeRange[1]);
              Double_t negative_percent = histos.at(ilcn->first)->Integral(b1Negative, b2Negative) / whole_integral;
              if(negative_percent > negativeCut){
                fChannelWords.at(ilcn->first) |= 1<<fCSSBits.GetBitIndex("qhs");
              }
              int b1Main = histos.at(ilcn->first)->FindBin(mainChargeRange[0]);
              int b2Main = histos.at(ilcn->first)->FindBin(mainChargeRange[1]);
              Double_t main_percent = histos.at(ilcn->first)->Integral(b1Main, b2Main) / whole_integral;
              if(main_percent < mainCut){
                fChannelWords.at(ilcn->first) |= 1<<fCSSBits.GetBitIndex("qhs");
              }
              int b1Second = histos.at(ilcn->first)->FindBin(secondChargeRange[0]);
              int b2Second = histos.at(ilcn->first)->FindBin(secondChargeRange[1]);
              Double_t second_percent = histos.at(ilcn->first)->Integral(b1Second, b2Second) / whole_integral;
              if(second_percent > secondCut){
                fChannelWords.at(ilcn->first) |= 1<<fCSSBits.GetBitIndex("qhs");
              }
              int b1Third = histos.at(ilcn->first)->FindBin(thirdChargeRange[0]);
              int b2Third = histos.at(ilcn->first)->FindBin(thirdChargeRange[1]);
              Double_t third_percent = histos.at(ilcn->first)->Integral(b1Third, b2Third) / whole_integral;
              if(third_percent > thirdCut){
                fChannelWords.at(ilcn->first) |= 1<<fCSSBits.GetBitIndex("qhs");
              }
              int b1Railed = histos.at(ilcn->first)->FindBin(railedChargeRange[0]);
              int b2Railed = histos.at(ilcn->first)->FindBin(railedChargeRange[1]);
              Double_t railed_percent = histos.at(ilcn->first)->Integral(b1Railed, b2Railed) / whole_integral;
              if(railed_percent > railedCut){
                fChannelWords.at(ilcn->first) |= 1<<fCSSBits.GetBitIndex("qhs");
              }
            }
        }

    void ChanSWStatusCalib::EndOfRun( DS::Run& run )
    {
        RAT::DBTable table( "CSS_RESULT" );
        table.SetI( "version", 2 );
        table.SetS( "index", "" );
        table.SetRunRange( fRunID, fRunID );
        table.SetPassNumber( -1 );
        table.SetZ( "production", false );
        table.SetS( "comment", dformat( "Created by ChanSWStatusCalib applied to %d", fRunID ) );
        // 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 );

        ApplyCHSCut();
        ApplyHighOccCut();
        ApplyECACut();
        ApplyPCACut();
        ApplyLowOccCut();
        // These lines refer to an unused test, they are commented out as the "qhs" flag is now being used elsewhere
        //ApplyChargeTimeCut( hQHS, "qhs" );
        //ApplyChargeTimeCut( hQHL, "qhl" );
        //ApplyChargeTimeCut( hQLX, "qlx" );
        //ApplyChargeTimeCut( hTAC, "tac" );
        ApplyQHSChargeCut(hQHS);

        table.SetIArray( "channelWord", fChannelWords );
        table.SaveAs( dformat("CSS_RESULT_%u.ratdb", run.GetRunID()) );

        // 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();
    }

    void
    ChanSWStatusCalib::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);
    }

}