#include <RAT/Processor.hh>
#include <RAT/DataQualityProc.hh>
#include <RAT/DQN16Proc.hh>
#include <RAT/DQTellieProc.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/PMTSet.hh>
#include <RAT/DS/UniversalTime.hh>
#include <RAT/DS/FitResult.hh>
#include <RAT/DS/FitVertex.hh>

#include <TH1I.h>
#include <TH1D.h>
#include <TF1.h>

#include <string>
#include <sstream>
#include <vector>
#include <bitset>

using namespace std;

namespace RAT{

    DQN16Proc::DQN16Proc() : DataQualityProc("dqn16proc"){
       fTH1INHitsTagged = NULL;
       fTH1INHitsUntagged = NULL;
       fTH1DTimeBetweenN16Events = NULL;
       fTH1DQHS= NULL;

       fN16Count = 0;
       fN16Events = 0;
    }

    DQN16Proc::~DQN16Proc(){
    }

    void DQN16Proc::BeginOfRun(DS::Run& run){
        DataQualityProc::BeginOfRun( run );
        //Load the database parameters for checks
        fDQChecks = DB::Get()->GetLink( "DQCHECKS", "n16" );
        //Get the DQChecks information
        //Get trigger bit
        fTrigCheckBit = fDQChecks->GetI("n16_trigger_bit");
        fCriteria["n16_trigger_bit"] = fTrigCheckBit;

        //Load NHit information
        fN16MeanNHitUpperLimit = fDQChecks->GetD("mean_nhits_upper");
        fCriteria["mean_nhits_upper"] = fN16MeanNHitUpperLimit;
        fN16MeanNHitLowerLimit = fDQChecks->GetD("mean_nhits_lower");
        fCriteria["mean_nhits_lower"] = fN16MeanNHitLowerLimit;
        fN16NHitSpreadUpperLimit = fDQChecks->GetD("nhits_spread_upper");
        fCriteria["nhits_spread_upper"] = fN16NHitSpreadUpperLimit;

        fN16FitTolerance =   fDQChecks->GetD("fit_position_max_diff");
        fCriteria["fit_position_max_diff"] = fN16FitTolerance;
        fMaxNumN16FitFails = fDQChecks->GetI("fit_position_max_fails");
        fCriteria["fit_position_max_fails"] = fMaxNumN16FitFails;


        //Load the run information from the database
        DBLinkPtr manipInfo = DB::Get()->GetLink("CALIB_COMMON_RUN_LEVEL","MANIP");
        vector<double> manipPos = manipInfo->GetDArray("position");
        fN16ManipPosition = TVector3(manipPos[0],manipPos[1],manipPos[2]);
        fAverageFitPosition  = TVector3(0.0,0.0,0.0);


        //Create the histograms
        stringstream titlestream;
        string title;
        titlestream << "_run_" << run.GetRunID();

        //NHits
        title = "N16 Trigged events Nhit"+titlestream.str();
        fTH1INHitsTagged = new TH1I(title.c_str(),title.c_str(),2000,0,2000);

        title = "Non N16 Trigged events Nhit"+titlestream.str();
        fTH1INHitsUntagged = new TH1I(title.c_str(),title.c_str(),2000,0,2000);
        //QHS
        title = "QHS_Histogram"+titlestream.str();
        fTH1DQHS= new TH1D(title.c_str(),title.c_str(),200,0,100);

        //Time between successive events
        title = "Time Between N16 events"+titlestream.str();
        fTH1DTimeBetweenN16Events = new TH1D(title.c_str(),title.c_str(),200,0,1000);

        //Initialize counters to 0
        fN16Count = 0;
        fFirstN16Event = true;
        fNumN16FitFails = 0;




    }

    void DQN16Proc::EndOfRun(DS::Run& run){
        //Write Histograms to ROOT file
        WriteToRoot(fTH1INHitsTagged);
        WriteToRoot(fTH1INHitsUntagged);
        WriteToRoot(fTH1DTimeBetweenN16Events);
        WriteToRoot(fTH1DQHS);
        info << "Wrote histograms to the root file"<<newline;
        info << "Number of N16 Events: "<<fN16Count<<newline;

        //Perform the gaussian fit
        if(fN16Count == 0){
            info << "Number of N16 Events is  0"<<newline;
            info << "Trigger Bit is: "<<fTrigCheckBit<<newline;
            fN16MeanNHit =  0;
            fN16NHitSpread = 9999;
        }
        else{
            fTH1INHitsTagged->Fit("gaus"); //Get the mean and sp
            TF1 * func = fTH1INHitsTagged->GetFunction("gaus");
            fN16MeanNHit = func->GetParameter(1);
            fN16NHitSpread = func->GetParameter(2);
            delete func;
            func = NULL;
        }
        //Checking mean NHit
        if(fN16MeanNHit < fN16MeanNHitLowerLimit || fN16MeanNHit > fN16MeanNHitUpperLimit){
            Update("n16MeanNHitCheck",false);
        }
        else{
            Update("n16MeanNHitCheck",true);
        }
        //Checking NHit spread
        if(fN16NHitSpread >  fN16NHitSpreadUpperLimit){
            Update("n16NHitSpreadCheck",false);
        }
        else{
            Update("n16NHitSpreadCheck",true);
        }
        //Checking event count
        if(fN16Count != fN16Events){
            Update("n16NumberOfEventsCheck",false);
        }
        else{
            Update("n16NumberOfEventsCheck",true);
        }
        if(fMaxNumN16FitFails > fNumN16FitFails){
            Update("n16FitPositionCheck",true);
        }
        else{
            Update("n16FitPositionCheck",false);
        }
        //Obtain the average Fit Position
        fAverageFitPosition.SetMag(fAverageFitPosition.Mag()/(fN16Count-fNumN16FitFails));
        stringstream outstream;
        outstream << "["<<fAverageFitPosition.X()<<", "<<fAverageFitPosition.Y()<<", "<<fAverageFitPosition.Z()<<"]";
        std::string averagePosOutput = outstream.str();
        AddToRATDB("mean_fit_position",averagePosOutput);
        //Add parameters to RATDB
        AddToRATDB("num_N16_events",fN16Count);
        AddToRATDB("mean_nhit",fN16MeanNHit);
        AddToRATDB("nhit_spread",fN16NHitSpread);
        AddToRATDB("num_events_failing_nhit_check",fNumN16FitFails);

        DataQualityProc::EndOfRun( run );
        //Delete the Histograms
        info << "Deleting the Histograms"<<newline;
        delete fTH1INHitsTagged;
        fTH1INHitsTagged = NULL;
        delete fTH1INHitsUntagged;
        fTH1INHitsUntagged = NULL;
        delete fTH1DTimeBetweenN16Events;
        fTH1DTimeBetweenN16Events = NULL;
        delete fTH1DQHS;
        fTH1DQHS= NULL;
    }

    Processor::Result DQN16Proc::DSEvent(DS::Run&, DS::Entry& ds)
    {
        for( size_t iEV = 0; iEV < ds.GetEVCount(); iEV++ )
            Event(ds, ds.GetEV(iEV));
        return OK;
    }

  Processor::Result DQN16Proc::Event( DS::Entry&, DS::EV& ev )
  {
    //Do the trigger cut
    bool bitCheck = ev.GetTrigType() & (1 <<fTrigCheckBit);
    //For events triggered by the N16 source
    if(bitCheck){
        fN16Count++;
        //Fill NHits Histogram
        fTH1INHitsTagged->Fill(ev.GetNhits());
        //Fill the QHS Histogram
        DS::CalPMTs evPMTs = ev.GetCalPMTs();
        for(unsigned int iPMT=0; iPMT<evPMTs.GetCount(); iPMT++){
            fTH1DQHS->Fill(evPMTs.GetPMT(iPMT).GetQHS());
        }

        //Check if first N16 event for calculating the time between events
        if(fFirstN16Event){
            fPrevEventTime = ev.GetUniversalTime();
            fFirstN16Event = false;
        }

        else{
            DS::UniversalTime deltaT = ev.GetUniversalTime()-fPrevEventTime;
            double timeInMilliS = ((deltaT.GetSeconds()*1e9)+deltaT.GetNanoSeconds())/1e6;
            fTH1DTimeBetweenN16Events->Fill(timeInMilliS);
        }

        //Now do fit results check
        try{
            const DS::FitResult& fitRes = ev.GetFitResult("N16Fit");
            int vertexCount = fitRes.GetVertexCount();
            //Should only be one vertex otherwise fail the check.
            if( vertexCount != 1){
                fNumN16FitFails++;
            }
            else{
                const DS::FitVertex& fitVert = fitRes.GetVertex(0);
                TVector3 fitPos =  fitVert.GetPosition();
                fAverageFitPosition = fAverageFitPosition + fitPos;
                if((fitPos-fN16ManipPosition).Mag() > fN16FitTolerance ){
                    fNumN16FitFails++;
                }
            }
        }
        catch(DS::DataNotFound){
            info << "No fit information for event: "<<ev.GetGTID()<<newline;
            fNumN16FitFails++;
        }

    }
    else{
        fTH1INHitsUntagged->Fill(ev.GetNhits());
    }
    return OK;
  }



}// namespace RAT