#include <TH1D.h>
#include <TH2D.h>
#include <TFile.h>
#include <TVector3.h>
#include <TMath.h>
#include <TGraph.h>
#include <TCanvas.h>

#include <RAT/Processor.hh>
#include <RAT/DataQualityProc.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DB.hh>
#include <RAT/DBLink.hh>
#include <RAT/DQPMTProc.hh>
#include <RAT/Log.hh>
#include <RAT/BitManip.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/ChanHWStatus.hh>

#include <vector>
#include <map>
#include <iterator>

using namespace std;

namespace RAT
{

  DQPMTProc::DQPMTProc() : DataQualityProc( "dqpmtproc" )
  {
  }

  DQPMTProc::~DQPMTProc()
  {
  }

  void DQPMTProc::BeginOfRun( DS::Run& run )
  {
    DataQualityProc::BeginOfRun( run );
    // Initialise channel type counts
    fNormalCount = 0;
    fOWLCount = 0;
    fLGCount = 0;
    fCalibCount = 0;
    fSpareCount = 0;
    fButtCount = 0;
    fNeckCount = 0;
    fFECDCount = 0;

    // Get variables from RATDB
    fGeoPsup = DB::Get()->GetLink( "NATIVE_GEO_DIMENSIONS", "natgeo_dimensions" );
    //fNeckRadius = 800.0;//fGeoPsup->GetD( "r_maxt" );
    //fPsupRadius = 8900.0;
    fNeckRadius = fGeoPsup->GetD("neck_inner_radius");
    fPsupRadius = fGeoPsup->GetD( "psup_radius" );

    fPanelInfo = DB::Get()->GetLink( "PANELINFO" );
    fPanelNumbers = fPanelInfo->GetIArray( "panel_number" );
    fPanelCount = fPanelNumbers.size();
    fPanelTypes = fPanelInfo->GetIArray( "panel_type" );

    // Get check values from RATDB, save as class members and to RATDB via JSON object
    fDQChecks = DB::Get()->GetLink( "DQCHECKS", "neutrino" );
    fGeneralCoverage = 0;
    fGeneralCovThresh = fDQChecks->GetD( "general_cov_thresh" );
    fCriteria[ "general_cov_thresh" ] = fGeneralCovThresh;
    fInCrateCovThresh = fDQChecks->GetD( "in_crate_cov_thresh" );
    fCriteria[ "in_crate_cov_thresh" ] = fInCrateCovThresh;
    fCrateCovThresh = fDQChecks->GetD( "crate_cov_thresh" );
    fCriteria[ "crate_cov_thresh" ] = fCrateCovThresh;
    fPanelCovThresh = fDQChecks->GetD( "panel_cov_thresh" );
    fCriteria[ "panel_cov_thresh" ] = fPanelCovThresh;

    fPMTInfo = DU::Utility::Get()->GetPMTInfo();
    fCHS = DU::Utility::Get()->GetChanHWStatus();
    fChannelCount = fPMTInfo.GetCount();

    fCrateCount = BitManip::GetCrate( fChannelCount );
    fCrateCoverage.resize( fCrateCount );

    fPMTCalStatuses.resize( fChannelCount );

    // Create Panel Map
    if( ( fPanelNumbers.size() == fPanelCount ) && ( fPanelTypes.size() == fPanelCount ) )
      {
        vector<int>::iterator iPanelNumber = fPanelNumbers.begin();
        vector<int>::iterator iPanelType = fPanelTypes.begin();
        while( iPanelNumber != fPanelNumbers.end() && iPanelType != fPanelTypes.end() )
          {
            int panelNumber = *iPanelNumber;

            fPanelProperties.type = *iPanelType;

            // Coverage threshold set at half the number of PMTs (rounding down if odd)
            if( fPanelProperties.type == 0 ) // Type 0 = S7
              fPanelProperties.coverageThresh = 3;
            if( fPanelProperties.type == 1 ) // Type 1 = S19
              fPanelProperties.coverageThresh = 9;
            if( fPanelProperties.type == 2 ) // Type 2 = T21
              fPanelProperties.coverageThresh = 10;
            if( fPanelProperties.type == 3 ) // Type 3 = T14
              fPanelProperties.coverageThresh = 7;
            if( fPanelProperties.type == 4 ) // Type 4 = T10
              fPanelProperties.coverageThresh = 5;

            fPanelProperties.coverageCount = 0;

            pair<map<int, PanelProp>::iterator, bool> addPanel;
            addPanel = fPanels.insert( pair<int, PanelProp>( panelNumber, fPanelProperties ) );
            if( addPanel.second == true)
              debug << "DQPMTProc::BeginOfRun: attempted to insert element that already exists in fPanels map" << newline;

            iPanelNumber++;
            iPanelType++;
          }
      }

    // Set PMT Type Counts
    for( vector<int>::iterator iLCN = fPMTCalStatuses.begin(); iLCN < fPMTCalStatuses.end(); iLCN++ )
      {
        int lcn = distance( fPMTCalStatuses.begin(), iLCN );
        if( fPMTInfo.GetType( lcn ) == DU::PMTInfo::NORMAL )
          fNormalCount++;
        else if( fPMTInfo.GetType( lcn ) == DU::PMTInfo::OWL )
          fOWLCount++;
        else if( fPMTInfo.GetType( lcn ) == DU::PMTInfo::LOWGAIN )
          fLGCount++;
        else if( fPMTInfo.GetType( lcn ) == DU::PMTInfo::CALIB )
          fCalibCount++;
        else if( fPMTInfo.GetType( lcn ) == DU::PMTInfo::SPARE )
          fSpareCount++;
        else if( fPMTInfo.GetType( lcn ) == DU::PMTInfo::BUTT )
          fButtCount++;
        else if( fPMTInfo.GetType( lcn ) == DU::PMTInfo::NECK )
          fNeckCount++;
        else
          debug << "DQPMTProc::BeginOfRun: PMT with lcn " << lcn << " did not mach a known PMT type" << newline;
      }

    // Calculate threshold area - size of neck
    // This code is preparation for a future feature
    const double PI = 3.14259265359;
    double capArea = GetCapArea( fPsupRadius, fNeckRadius );
    fNormalPMTsArea = ( 4*PI*fPsupRadius*fPsupRadius ) - capArea;
    fCapCount = static_cast<double>( fNormalCount )/fNormalPMTsArea*capArea;
    fNeckAreaCovThresh = fCapCount/static_cast<double>( fNormalCount );

    // Book histograms/TGraphs
    fCrateOccupancyMap = new TH2D( "TH2D_CrateCoverageMap",
                                    "PMT coverage map after PMT Calibrations",
                                    320, 0.0, 20.0, 32, 0.0, 32.0 );
    fCrateOccupancyMap->SetDirectory(0);
    fCrateOccupancyMap->SetStats(0);
    fUnCalCrateOccupancyMap = new TH2D( "TH2D_UnCalCrateCoverageMap",
                                         "PMT coverage map before PMT Calibrations",
                                         320, 0.0, 20.0, 32, 0.0, 32.0 );
    fUnCalCrateOccupancyMap->SetDirectory(0);
    fUnCalCrateOccupancyMap->SetStats(0);

    fUncalPMTs = new TH1D( "TH1D_UnCalCrateCoverageMap", "PMT coverage map before PMT Calibrations", 10000, 0, 10000);
    fUncalPMTs->SetDirectory(0);
    fUncalPMTs->SetStats(0);

    fcalPMTs = new TH1D( "TH1D_CalCrateCoverageMap", "PMT coverage map after PMT Calibrations", 10000, 0, 10000);
    fcalPMTs->SetDirectory(0);
    fcalPMTs->SetStats(0);
  }

  double DQPMTProc::GetCapArea( const double r, const double a ) const
  {
    const double PI = 3.14259265359;
    double result = (-1) * TMath::Sqrt( r*r - a*a );
    result += r;
    result *= 2*PI*r;
    return result;
  }

  void DQPMTProc::FillCrateCoverageMap( const int lcn, TH2D* h1 ) const
  {
    // Fill crate view histogram
    int crate = BitManip::GetCrate( lcn );
    int card = BitManip::GetCard( lcn );
    int channel = BitManip::GetChannel( lcn );
    double x = crate + static_cast<double>( card )/16.0;
    h1->Fill( x, channel );
  }

  void DQPMTProc::FillCrateCoverageMap( const int lcn, TH1D* h1) const
  {
    h1->SetBinContent(lcn, 1.0);
  }

  void DQPMTProc::UpdateCoverage( const int lcn, const DU::PMTInfo& pmtInfo)
  {
      if(fCHS.IsChannelOnline(lcn) && fCHS.IsTubeOnline(lcn) && fCHS.IsDAQEnabled(lcn)){
          fGeneralCoverage++;

          int crate = BitManip::GetCrate( lcn );
          fCrateCoverage[ crate ]++;

          int panelNumber = pmtInfo.GetPanelNumber( lcn );
          map<int, PanelProp>::iterator iPanel = fPanels.find( panelNumber );
          if( panelNumber != iPanel->first )
          {
              debug << "DQPMTProc::UpdateCoverage: panel number from iterator (" << panelNumber;
              debug << ") does not equal panel number from DU::PMTInfo::GetPanelNumber( lcn )" << newline;
          }
          else
              iPanel->second.coverageCount++;
      }
  }

  void DQPMTProc::RunCoverageChecks()
  {
    // Check 1: General detector coverage
    double detectorCoverage = fGeneralCoverage/static_cast<double>( fNormalCount )*100.0; // normalised by normal count
    AddToRATDB("overall_detector_coverage",detectorCoverage);
    if( detectorCoverage > fGeneralCovThresh )
      {
        detail << "DQPMTProc::RunCoverageChecks: run " << GetRunID() << " passed general detector coverage check" << newline;
        fCheckResult = true;
      }
    else
      {
        detail << "DQPMTProc::RunCoverageChecks: run " << GetRunID() << " failed general detector coverage check" << newline;
        detail << " --> general coverage (" << detectorCoverage << "%) was below the " << fGeneralCovThresh << "% general coverage threshold " << newline;
      }
    Update( "general_coverage", fCheckResult );

    // Check 2: Crate coverage check
    int failedCrates=0;
    for( vector<double>::iterator iCrateCoverage = fCrateCoverage.begin(); iCrateCoverage != fCrateCoverage.end(); iCrateCoverage++)
      {
        *iCrateCoverage /= 512.0; // Normalised by number of channels per crate
        *iCrateCoverage *= 100.0; // Turn into %
        if( *iCrateCoverage < fInCrateCovThresh )
          {
            detail << "DQPMTProc::RunCoverageChecks: crate " << distance( fCrateCoverage.begin(), iCrateCoverage ) << " has coverage ";
            detail << *iCrateCoverage << "%, which is below the " << fInCrateCovThresh << "% threshold" << newline;
            failedCrates++;
          }
      }
    AddToRATDB("crates_failing_coverage",failedCrates);
    double crateCoverage = 100.0 * ( 1.0 - ( static_cast<double>( failedCrates) / static_cast<double>( fCrateCount ) ) );
    AddToRATDB("crates_coverage_percentage",crateCoverage);
    if ( crateCoverage >= fCrateCovThresh )
      {
        detail << "DQPMTProc::RunCoverageChecks: run " << GetRunID() << " passed crate coverage check" << newline;
        fCheckResult = true;
      }
    else
      {
        detail << "DQPMTProc::RunCoverageChecks: run " << GetRunID() << " failed crate coverage check" << newline;
        detail << " --> " << failedCrates << " crates had coverage below the " << fInCrateCovThresh << "% in-crate coverage threshold " << newline;
      }
    Update( "crate_coverage", fCheckResult );

    // Check 3: Panel coverage check
    int failedPanels = 0;
    int checkPMTCount = 0;
    for( map<int, PanelProp>::iterator iPanel = fPanels.begin(); iPanel != fPanels.end(); iPanel++ )
      {
        int coverage = iPanel->second.coverageCount;
        int threshold = iPanel->second.coverageThresh;
        checkPMTCount += coverage;
        if ( coverage <= threshold )
          {
            detail << "DQPMTProc::RunCoverageChecks: panel " << iPanel->first << " has coverage ";
            detail << coverage << " PMTs, which is below the " << threshold << " PMT threshold" << newline;
            failedPanels++;
          }
      }
    if( checkPMTCount != fGeneralCoverage )
      {
        debug << "DQPMTProc::RunCoverageChecks: panel check counted " << checkPMTCount << " panels which does not match ";
        debug << fGeneralCoverage << " panels counted by General Detector Coverage check" << newline;
      }
    AddToRATDB("number_of_panels_failing_coverage",failedPanels);
    double panelCoverage = 100.0 * ( 1.0 - ( static_cast<double>( failedPanels ) / static_cast<double>( fPanelCount ) ) );
    AddToRATDB("percentage_of_panels_passing_coverage",panelCoverage);
    if( panelCoverage >= fPanelCovThresh )
      {
        detail << "DQPMTProc::RunCoverageChecks: run " << GetRunID() << " passed panel coverage check" << newline;
        fCheckResult = true;
      }
    else
      {
        detail << "DQPMTProc::RunCoverageChecks: run " << GetRunID() << " failed panel coverage check" << newline;
        detail << " --> " << failedPanels << " panels had coverage below their panel coverage threshold" << newline;
      }
    Update( "panel_coverage", fCheckResult );

    // Check 4: Neck sized area, no coverage check -> to be added
  }

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

  Processor::Result DQPMTProc::Event( DS::Entry&, DS::EV& ev )
  {
    fEventCount++;
    int eventID = ev.GetGTID();

    const size_t PMTAllUnCalCount = ev.GetUncalPMTs().GetAllCount();
    if( PMTAllUnCalCount == 0 )
      debug << "DQPMTProc::Event: no un-calibrated PMTs were found for event " << eventID << newline;
    size_t hitPMTCount = 0;
    const size_t PMTAllCalCount = ev.GetCalPMTs().GetAllCount();
    if( PMTAllCalCount == 0 )
      debug << "DQPMTProc::Event: No calibrated PMTs were found for event " << eventID << newline;
    size_t calibratedPMTCount = 0;

    vector<int>::iterator iPMTCalStatus = fPMTCalStatuses.begin();
    while( iPMTCalStatus != fPMTCalStatuses.end() )
      {
        for( size_t ipmt = 0; ipmt < PMTAllUnCalCount; ipmt++ )
          {
            // Has the PMT been hit? i.e. does a DS::GetPMTALLUnCal* object exist
            try
              {
                const DS::PMTUncal& pmtuncal = ev.GetUncalPMTs().GetAllPMT( ipmt );
                unsigned int lcn = pmtuncal.GetID();
                int moveIterator = lcn - distance( fPMTCalStatuses.begin(), iPMTCalStatus );
                iPMTCalStatus += moveIterator; // *iPMTCalStatus now point to correct entry
                if( lcn != distance( fPMTCalStatuses.begin(), iPMTCalStatus ) )
                  {
                    debug << "DQPMTProc::Event: lcn from iterator (" << distance( fPMTCalStatuses.begin(), iPMTCalStatus );
                    debug << ") does not equal lcn from DS::UncalPMT::GetID() (" << lcn << ")" << newline;
                  }

                if( *iPMTCalStatus == 0 ) // PMT has not yet been hit
                  {
                    *iPMTCalStatus = 1 ; // Change to status 1 (means PMT has been hit at least once)
                    debug << "DQPMTProc::Event: changing PMT " << lcn << " calibration status to 1" << newline;
                  }

                hitPMTCount++;
                FillCrateCoverageMap( lcn, fUnCalCrateOccupancyMap );
                FillCrateCoverageMap( lcn, fUncalPMTs );

                bool isCalibrated = false;
                // Has PMT been calibrated? i.e. does a DS::GetPMTALLCal* object exist
                try
                  {
                    const DS::PMTCal& pmtcal = ev.GetCalPMTs().GetAllPMT( ipmt );
                    isCalibrated = true;
                    if( pmtcal.GetID() == lcn ) // Check for the initial hit pmt, if it was also calibrated
                      {
                        if( *iPMTCalStatus == 1 ) // PMT has been hit
                          {
                            *iPMTCalStatus = 3; // Change to status 3 (means each time PMT has been hit it has been calibrated)
                            debug <<  "DQPMTProc::Event: changing PMT " << lcn << " calibration status to 3" << newline;
                          }
                        calibratedPMTCount++;
                        FillCrateCoverageMap( lcn, fCrateOccupancyMap );
                        FillCrateCoverageMap( lcn, fcalPMTs );
                      }
                  }
                catch (const out_of_range& oor)
                  {
                    debug << "DQPMTProc::Event: No object of type DS::GetPMTALLCal* for PMT " << ipmt << newline;
                  }

                if( !isCalibrated && *iPMTCalStatus == 3 ) // PMT Has been calibrated previously but has not been calibrated for this event
                  {
                    *iPMTCalStatus = 2; // Change to status 2 (means PMT has not been calibrated every time it was hit)
                    debug << "DQPMTProc::Event: changing PMT " << lcn << " calibration status to 2" << newline;
                  }
              }
            catch (const out_of_range& oor)
              {
                debug << "DQPMTProc::Event: no object of type DS::GetPMTAllUnCal for pmt " << ipmt << newline;
              }
          }
        iPMTCalStatus = fPMTCalStatuses.end();
      }

    if( hitPMTCount != PMTAllUnCalCount )
      {
        debug << "DQPMTProc::Event: calculated number of hit PMTs (" << hitPMTCount;
        debug << ") differs from PMTAllUnCalCount (" << PMTAllUnCalCount << ") for this event" << newline;
      }
    if( calibratedPMTCount != PMTAllCalCount )
      {
        debug << "DQPMTProc::Event: calculated number of calibrated PMTs (" << calibratedPMTCount;
        debug << ") differs from PMTAllCalCount (" << PMTAllCalCount << ") for this event" << newline;
      }

    return OK;
}

  void DQPMTProc::EndOfRun( DS::Run& run )
  {
    // Loop over PMTs, add records to histograms/scatter plots for PMTs that pass calibration criteria
    int ipmt = 0;
    for( vector<int>::iterator iPMTCalStatus = fPMTCalStatuses.begin(); iPMTCalStatus !=fPMTCalStatuses.end(); iPMTCalStatus++ )
      {
        int lcn = distance( fPMTCalStatuses.begin(), iPMTCalStatus );
        if( *iPMTCalStatus == 3 && fPMTInfo.GetType( lcn ) == DU::PMTInfo::NORMAL ) // PMT passes calibration criteria and is Normal
          {
            UpdateCoverage( lcn, fPMTInfo);
            TVector3 pos = fPMTInfo.GetPosition( lcn );
            double phi = pos.Phi();
            double theta = pos.Theta();
            fPMTPhi.push_back( phi );
            fPMTTheta.push_back( theta );
            ipmt++;
          }
      }
    RunCoverageChecks();

    // write hist to root file
    fGeoCoverageMap = new TGraph( fPMTPhi.size(),
                                     &fPMTPhi[0],
                                     &fPMTTheta[0] );
    fGeoCoverageMap->SetName( "TGraph_GeoCoverageMap" );

    if( fOutputImages ) // Only create image files if this has been set to true
      {
        Plot2DHist( fCrateOccupancyMap, "Crate", "Channel" );
        Plot2DHist( fUnCalCrateOccupancyMap, "Crate", "Channel" );
        Plot2DScatter( fGeoCoverageMap, "Phi", "Theta", 21, 38, 0.35 );
      }

    WriteToRoot( fCrateOccupancyMap, "fCrateOccupancyMap" );
    WriteToRoot( fUnCalCrateOccupancyMap, "fUnCalCrateOccupancyMap" );
    WriteToRoot( fGeoCoverageMap, "fGeoCoverageMap" );

    TH2D pmtCrateCoverageMap = CrateMapFromChannelHisto(*fcalPMTs); //create the 2D hist
    TH2D pmtUncalCrateCoverageMap = CrateMapFromChannelHisto(*fUncalPMTs);
    TCanvas *C1 = new TCanvas("C1");
    string filename = pmtCrateCoverageMap.GetName();
    //fCrateOccupancyMap->SetXTitle("Crate");
    //fCrateOccupancyMap->SetYTitle("Channel");
    filename = filename+".png";
    C1->cd();
    //fCrateOccupancyMap->Draw();
    pmtCrateCoverageMap.Draw("colz");
    C1->Print( filename.c_str() );

    TCanvas *C2 = new TCanvas("C2");
    filename = pmtUncalCrateCoverageMap.GetName();
    filename = filename+".png";
    //fUnCalCrateOccupancyMap->SetXTitle("Crate");
    //fUnCalCrateOccupancyMap->SetYTitle("Channel");
    //fUnCalCrateOccupancyMap->Draw();
    pmtUncalCrateCoverageMap.Draw("colz");
    C2->Print( filename.c_str() );

    TCanvas *C3 = new TCanvas("C3"); //new canvas for TGraph fGeoCoverageMap
    filename = fGeoCoverageMap->GetName();
    fGeoCoverageMap->SetTitle("Geo Coverage Map; Phi; Theta");
    filename = filename+".png";
    C3->cd();
    fGeoCoverageMap->Draw("ap");
    C3->Print( filename.c_str() );
    delete C1;
    C1 = NULL;
    delete C2;
    C2 = NULL;
    delete C3;
    C3 = NULL;

    delete fUncalPMTs;
    fUncalPMTs = NULL;
    delete fcalPMTs;
    fcalPMTs = NULL;
    delete fCrateOccupancyMap;
    fCrateOccupancyMap = NULL;
    delete fUnCalCrateOccupancyMap;
    fUnCalCrateOccupancyMap = NULL;
    delete fGeoCoverageMap;
    fGeoCoverageMap = NULL;

    DataQualityProc::EndOfRun( run );
  }

  TH2D DQPMTProc::CrateMapFromChannelHisto(TH1D channelFails){
    char histoName [100];
    sprintf(histoName,"%s",channelFails.GetTitle());
    TH2D outputHisto(histoName,histoName,280,0,1,100,0,1);
    int crateGap = 10;
    int xOffset = 10;
    int yOffset = 10;

    for(unsigned int i=0; i<fPMTInfo.GetCount(); i++){
      int crateNum = crateNumber(i);
      int cardNum = cardNumber(i);
      int channelNum = channelNumber(i);
      int ybin = channelNum+yOffset;
      int xbin = cardNum+(crateNum*(16+crateGap))+xOffset;
      if( crateNum >= 9){
        ybin += 50;
        xbin = cardNum+((crateNum-9)*(16+crateGap))+xOffset;
      }
      //16 cards per crate
      int bin = outputHisto.GetBin(xbin,ybin);
      outputHisto.SetBinContent(bin,channelFails.GetBinContent(i)*10000+1);
    }
    outputHisto.SetStats(false);
    return outputHisto;
  }

  int DQPMTProc::crateNumber(int lcn){
    return lcn/512;
  }

  int DQPMTProc::cardNumber(int lcn){
    return ((lcn-(512*crateNumber(lcn)))/32);
  }

  int  DQPMTProc::channelNumber(int lcn){
    return lcn-(32*cardNumber(lcn))-(512*crateNumber(lcn));
    }
}// namespace RAT