// G4 stuff
#include <G4UImessenger.hh>
#include <G4PhysicalConstants.hh>

// Root stuff
#include <TFile.h>
#include <TVector.h>
#include <TMath.h>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TH1F.h>
#include <TGraph2D.h>
#include <TImage.h>
#include <TView3D.h>
#include <TCanvas.h>
#include <TStyle.h>

// RAT stuff
#include <RAT/PCATellieMonitorProc.hh>
#include <RAT/PCAHists.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/LightPathCalculator.hh>
#include <RAT/DU/ChanHWStatus.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/FlatMap.hh>
#include <RAT/PhysicsUtil.hh>

// C++ stuff
#include <sys/stat.h>
#include <math.h>
#include <stdio.h>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <limits.h>

namespace RAT {

  PCATellieMonitorProc::PCATellieMonitorProc() : Processor("PCAMonitorProc")
  {
    // Set the total number of PMTs = nCrates * nCards * nChannels = 19*16*32
    fTotPMTs = 19 * 16 * 32;
    // Set up the variables that keep track of the minimum and maximum NHit
    // over all events
    // The maximum number of hits per event is the total number of PMTs, so
    // fNHitsMin will be smaller than this. Initialize it to the maximum + 1.
    fNHitsMin = fTotPMTs + 1;
    // The minimum number of hits per events is 0, therefore
    // fNHitsMax will be larger than this. Initialize it to -1
    fNHitsMax = -1;
    fTransitTimes.resize(fTotPMTs);
    fLEDIniDir.resize(fTotPMTs);
    fHitCountCentral.resize(fTotPMTs);
    fHitCount.resize(fTotPMTs);
    fEventCounter = 0;
    for (int i = 0; i < fTotPMTs; i++) {
      fTransitTimes[i] = -9999.0;
      fHitCountCentral[i] = 0;
      fHitCount[i] = 0;
    }
    // Set the flag that indicates that the transittimes were calculated
    fRunCalcTransitTimes = FALSE;
    // Set the flag that indicates that run-specific checks have been performed
    fRunCheckOK = FALSE;
    // Grab the default PCA settings from the PCA_GENERATION table
    DBLinkPtr PCAbank = DB::Get()->GetLink("PCA_GENERATION");
    fSelectAngle = PCAbank->GetD("select_angle");
    fTrigType = PCAbank->GetI("trig_type");
    // Set up the histograms
    // Axis range is large, will be adjusted when the histograms are saved
    fQHSCal = new TH1F("QHScal", "QHSCal", 4000, -1000, 3000);
    fQHLCal = new TH1F("QHLcal", "QHLCal", 4000, -1000, 3000);
    fQHSUnCal = new TH1F("QHSuncal", "QHSUnCal", 4095, 0, 4095);
    fQHLUnCal = new TH1F("QHLuncal", "QHLUnCal", 4095, 0, 4095);
    fTimeCal = new TH1F("Timecal", "TimeCal", 4000, -1000, 3000);
    fTimeUnCal = new TH1F("Timeuncal", "TimeUnCal", 4095, 0, 4095);
    fTimeResidual = new TH1F("Timeresidual", "TimeResidual", 4000, -1000, 3000);
  }

  void PCATellieMonitorProc::BeginOfRun(DS::Run&)
  {
    // Reset the flag that indicates that the transittimes where calculated
    fRunCalcTransitTimes = FALSE;
    // Reset the flag that indicates that run-specific checks have
    // been performed
    fRunCheckOK = FALSE;
    // Grab the LED position (fInputFiber) from the data-base
    DBLinkPtr ledInfo;
    fLEDPositionX = -9999.0;
    fLEDPositionY = -9999.0;
    fLEDPositionZ = -9999.0;
    fLEDDirX = -9999.0;
    fLEDDirY = -9999.0;
    fLEDDirZ = -9999.0;
    // Try to retrieve the position and direction of the fiber from the database
    std::ostringstream oss;
    oss.str("");
    // Take care of zero-padding
    oss << "FT" << std::setw(3) << std::setfill('0') << fInputFiber;
    // hard-coded to position 'A' at this time.
    // FIXME: in future, we need to find out which fibers are used
    // from an 'installed/used' entry in the database, not from user-input
    oss << "A";
    try {
      ledInfo =  DB::Get()->GetLink("FIBRE", oss.str());
      fLEDPositionX = ledInfo->GetD("x");
      fLEDPositionY = ledInfo->GetD("y");
      fLEDPositionZ = ledInfo->GetD("z");
      fLEDDirX = ledInfo->GetD("u");
      fLEDDirY = ledInfo->GetD("v");
      fLEDDirZ = ledInfo->GetD("w");
    }
    catch (DBNotFoundError &e) {
      info << "PCATellieMonitorProc(): Fiber " << oss.str()
           << " does not exist." << endl;
      // We should not proceed if we were unable to retrieve the fiberposition
      // Pass on the exception
      throw &e;
    }
  }

  void PCATellieMonitorProc::SetI(const std::string& param, const int value)
  {
    if (param == std::string("fiber")) {
      if (value >= 0 && value <= 115)
        fInputFiber = value;
      else
        throw ParamInvalid(param, "value must be between 0 and 115");
    } else {
      throw ParamUnknown(param);
    }
  }

  PCATellieMonitorProc::~PCATellieMonitorProc()
  {
    // Save the accumulated time and charge histograms
    PCATellieMonitorProc::Save1DHisto(fQHSCal, "QHS_cal [cap]");
    PCATellieMonitorProc::Save1DHisto(fQHLCal, "QHL_cal [cap]");
    PCATellieMonitorProc::Save1DHisto(fQHSUnCal, "QHS_uncal [counts]");
    PCATellieMonitorProc::Save1DHisto(fQHLUnCal, "QHL_uncal [counts]");
    PCATellieMonitorProc::Save1DHisto(fTimeCal, "Time_ecacal [ns]");
    PCATellieMonitorProc::Save1DHisto(fTimeUnCal, "Time_uncal [counts]");
    PCATellieMonitorProc::Save1DHisto(fTimeResidual, "Time_residual [ns]");
    // The named errors in a map
    std::map<std::string, bool> flags;
    flags["Overall failed PCA TELLIE run"] = FALSE;
    flags["Number of events too low"] = FALSE;
    flags["Average number of hits per PMT is too low"] = FALSE;
    flags["Average occupancy per PMT is too high"] = FALSE;
    flags["The spread of NHit/event is too large"] = FALSE;
    flags["The NHit/event is not constant"] = FALSE;
    // Set up some counters
    int pmtCounter = 0;
    int pmtSum = 0;
    // Set up the variables that will hold the maximum and minimum
    // number of hits per PMT.
    int maxHits = -9000;
    int minHits = 9000;
    // More histograms!
    TH1F *hitsPerPMT = new TH1F("HitsPerPMT", "HitsPerPMT", 200000, 0, 200000);
    TH1F *hitsPerPMTAll = new TH1F("HitsPerPMTAll",
                                   "HitsPerPMTAll", 100000, 0, 100000);
    TH1F *occPerPMT = new TH1F("OccPerPMT", "OccPerPMT", 50000, 0, 50);
    // TGraph2D for the detector flatmaps
    TGraph2D *coverage = new TGraph2D();
    TGraph2D *coverageSelect = new TGraph2D();
    // We need PMTinfor and hardware status
    const DU::PMTInfo& pmtList = DU::Utility::Get()->GetPMTInfo();
    const DU::ChanHWStatus& chs = DU::Utility::Get()->GetChanHWStatus();
    // Set the LED position
    TVector3 ledPos;
    ledPos.SetX(fLEDPositionX);
    ledPos.SetY(fLEDPositionY);
    ledPos.SetZ(fLEDPositionZ);
    // Time to loop over all channels, fill the flatmaps and do some
    // bookkeeping
    int graphIndex = 0;
    for (int i = 0; i < fTotPMTs; i++) {
      TVector3 pmtPos = pmtList.GetPosition(i);
      TVector2 pmtPosFlat;
      // Use the FlatMap function to get the 2D position
      RAT::SphereToIcosahedron(pmtPos, pmtPosFlat, 1, 2.12);
      // We are not plotting the offline channels
      // Also check for NAN positions
      if (chs.IsTubeOnline(i) && (pmtPosFlat.X() == pmtPosFlat.X())) {
        coverage->SetPoint(graphIndex, pmtPosFlat.X(), pmtPosFlat.Y(),
                           fHitCount[i]);
        coverageSelect->SetPoint(graphIndex, pmtPosFlat.X(),
                                 pmtPosFlat.Y(), fHitCountCentral[i]);
        graphIndex++;
      } else {
        continue;
      }
      hitsPerPMTAll->Fill(fHitCount[i]);
      if (fHitCountCentral[i] > 0) {
        if (fHitCountCentral[i] > maxHits) maxHits = fHitCountCentral[i];
        if (fHitCountCentral[i] < minHits) minHits = fHitCountCentral[i];
        hitsPerPMT->Fill(fHitCountCentral[i]);
        pmtCounter++;
        pmtSum = pmtSum + fHitCountCentral[i];
        occPerPMT->Fill(fHitCountCentral[i]
                        / static_cast<double>(fEventCounter));
      }
    }
    // Some more calculations are needed
    double averageOccupancy = pmtSum
      / (static_cast<double>(pmtCounter)
         * fEventCounter);
    double averageNumberOfHits = pmtSum / (static_cast<double>(pmtCounter));
    // Save the histograms
    PCATellieMonitorProc::Save1DHisto(hitsPerPMT, "Number of hits");
    PCATellieMonitorProc::Save1DHisto(hitsPerPMTAll, "Number of hits");
    PCATellieMonitorProc::Save1DHisto(occPerPMT, "PMT occupancy");
    // Save the flatmaps
    PCATellieMonitorProc::SaveTGraph2D(coverage, "coverage");
    PCATellieMonitorProc::SaveTGraph2D(coverageSelect, "coverageselect");
    // Set up Timage and TCanvas for remaining, special plots
    TImage *img = TImage::Create();
    TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);
    char imagename[30];
    // Time to make some decisions on the quality of this run!
    // 1) Are there enough events?
    // To start with, we would like a minimum of 10000 events
    if (fEventCounter < 10000) {
      // Not enough events!
      flags["Overall failed PCA TELLIE run"] = TRUE;
      flags["Number of events too low"] = TRUE;
    }
    // 2) Are there enough events per PMT on average?
    // For a good PCA, the central region PMTs should have ~500 hits
    if (averageNumberOfHits < 500) {
      // Not enough hits!
      flags["Overall failed PCA TELLIE run"] = TRUE;
      flags["Average number of hits per PMT is too low"] = TRUE;
    }
    // 3) Is the occupancy too high on average?
    // For good pca, the occupancy should be less than 5%
    // in central region
    if (averageOccupancy > 0.05) {
      // Occupancy is too high!
      flags["Overall failed PCA TELLIE run"] = TRUE;
      flags["Average occupancy per PMT is too high"] = TRUE;
    }
    // 4) Did the nhit change too much from event to event?
    // This would be erratic behaviour, test by looking at width of
    // nhit distribution.
    TH1F *nHits = new TH1F("nHits", "nHits", fTotPMTs, 0 , fTotPMTs);
    for (size_t i = 0; i < fNHitCheck.size(); i++) {
      nHits->Fill(fNHitCheck[i]);
    }
    TF1 * gaussian;
    double max = nHits->GetBinCenter(nHits->GetMaximumBin());
    gaussian = new TF1("gaussian", "gaus", max - 50, max + 50);
    nHits->Fit(gaussian, "Rq");
    double nhitSpread = gaussian->GetParameter(2);
    double avNhit = gaussian->GetParameter(1);
    if (nhitSpread > 25.0) {
      // NHit spread is too high
      flags["Overall failed PCA TELLIE run"] = TRUE;
      flags["The spread of NHit/event is too large"] = TRUE;
    }
    c1->cd();
    c1->Clear();
    nHits->SetAxisRange(nHits->FindFirstBinAbove(0),
                        nHits->FindLastBinAbove(0));
    nHits->GetXaxis()->SetLabelSize(0.03);
    nHits->GetYaxis()->SetLabelSize(0.03);
    nHits->GetXaxis()->SetTitle("Nhit");
    nHits->GetYaxis()->SetTitle("Counts");
    nHits->SetTitle("");
    nHits->SetStats(0);
    nHits->Draw();
    gaussian->Draw("SAME");
    img->Clear();
    img->FromPad(c1);
    snprintf(imagename, sizeof(imagename), "NHits-%i.png", fRunID);
    img->WriteImage(imagename);
    // 5) Did the nhit slip?
    // This could be because of the heating up of the board.
    // This is only an issue for PCA if it means that some events
    // are above maximum occupancy...
    // Need to fill these arrays to feed to TGraph
    size_t numPoints = fNHitCheck.size();
    double *xarray = new double[numPoints];
    double *yarray = new double[numPoints];
    for (size_t i = 0; i < numPoints; i++) {
      xarray[i] = i;
      yarray[i] = fNHitCheck[i];
    }
    TGraphErrors *nhitGraph = new TGraphErrors(numPoints, xarray, yarray);
    TF1 *polynomial = new TF1("pol1", "pol1", 0, numPoints);
    nhitGraph->Fit(polynomial, "Rq");
    if (TMath::Abs(polynomial->GetParameter(1))
        - (polynomial->GetParError(1) / 3.0) < 0) {
      // NHit is not stable
      flags["Overall failed PCA TELLIE run"] = TRUE;
      flags["The NHit/event is not constant"] = TRUE;
    }
    c1->cd();
    c1->Clear();
    c1->SetLogy(0);
    gStyle->SetOptStat(0);
    nhitGraph->GetXaxis()->SetTitle("Event number");
    nhitGraph->GetYaxis()->SetTitle("NHit");
    nhitGraph->SetTitle("");
    nhitGraph->Draw();
    img->Clear();
    img->FromPad(c1);
    snprintf(imagename, sizeof(imagename), "NHitvsEvent-%i.png", fRunID);
    img->WriteImage(imagename);
    // Clean up
    if (c1) c1->Close();
    img->Delete();
    // Now write this out to a file that can be read in by the python
    // analysis script.
    RAT::DBTable PCAMonitorTable("pca_monitor");
    PCAMonitorTable.SetD("nhit_spread", nhitSpread);
    PCAMonitorTable.SetD("av_occ", averageOccupancy);
    PCAMonitorTable.SetD("av_hit", averageNumberOfHits);
    PCAMonitorTable.SetD("av_nhit", avNhit);
    PCAMonitorTable.SetD("nevents", fEventCounter);
    PCAMonitorTable.SetI("flag_overall_fail",
                         flags["Overall failed PCA TELLIE run"]);
    PCAMonitorTable.SetI("fiber_number", fInputFiber);
    PCAMonitorTable.SetI("flag_nevents",
                         flags["Number of events too low"]);
    PCAMonitorTable.SetI("flag_av_hit",
                         flags["Average number of hits per PMT is too low"]);
    PCAMonitorTable.SetI("flag_av_occ",
                         flags["Average occupancy per PMT is too high"]);
    PCAMonitorTable.SetI("flag_nhit_spread",
                         flags["The spread of NHit/event is too large"]);
    PCAMonitorTable.SetI("flag_nhit_walk",
                         flags["The NHit/event is not constant"]);
    char filename[60];
    snprintf(filename, sizeof(filename), "PCAMonitorStatus-%i.ratdb", fRunID);
    PCAMonitorTable.SaveAs(filename);

    // Print a summary to the screen.
    std::cout << "" << std::endl;
    std::cout << "--------------" << std::endl;
    std::cout << "--- PCATellieMonitorProc(): SUMMARY " << std::endl;
    if (flags["Overall failed PCA TELLIE run"] == TRUE) {
      std::cout << "--- Run quality for run " << fRunID
                << " is BAD: " << std::endl;
      std::cout << "--- Details: " << std::endl;
      std::map<std::string, bool>::iterator it = flags.begin();
      while (it != flags.end()) {
        std::cout << "--- " << it->first << " :: "
                  << it->second << std::endl;
        it++;
      }
    } else {
      std::cout << "--- Run quality for run " << fRunID
                << " is GOOD" << std::endl;
    }
    std::cout << "--------------" << std::endl;
    std::cout << "" << std::endl;
  }

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

  Processor::Result PCATellieMonitorProc::Event(DS::Entry& ds, DS::EV& ev)
  {
    if (ev.GetTrigType() != fTrigType) return OK;   // Require EXTASYNC trigger
    fEventCounter++;   // This is an ASYNC event, therefore: count it!
    // Set the LED position
    TVector3 ledPos;
    TVector3 ledDir;
    ledPos.SetX(fLEDPositionX);
    ledPos.SetY(fLEDPositionY);
    ledPos.SetZ(fLEDPositionZ);
    ledDir.SetX(fLEDDirX);
    ledDir.SetY(fLEDDirY);
    ledDir.SetZ(fLEDDirZ);
    TVector3 injectedPos;
    TVector3 pmtPos;
    // Link to PMT information (position, type, etc)
    const DU::PMTInfo& pmtList = DU::Utility::Get()->GetPMTInfo();
    // lightpath instance is needed to retrieve the transittime
    // for the refracted path
    DU::LightPathCalculator lightPath = DU::Utility::Get()
      ->GetLightPathCalculator();
    // Use the following if-section to perform operarions that need to be
    // done only once per run and cannot be done in BeginOfRun
    if (fRunCheckOK !=TRUE) {
      // No additional checks here for now
      fRunID = ds.GetRunID();
      fRunCheckOK = TRUE;
    }
    // Calculating the transittime from the source position to the PMTs is slow,
    // therefore, do it at the start of the run for all PMTs.
    // We are assuming that there is one fiber firing per run.
    if (fRunCalcTransitTimes != TRUE) {
      for (int i = 0; i < fTotPMTs; i++) {
        // First, set the PMT position
        pmtPos = pmtList.GetPosition(i);
        // FIXME: hardcoded wavelength for now.
        double lambda = 505.0;
        lightPath.CalcByPosition(ledPos, pmtPos,
                                RAT::util::WavelengthToEnergy(lambda
                                                              * CLHEP::nm),
                                30.0);
        double distInInnerAV = lightPath.GetDistInInnerAV();
        double distInAV = lightPath.GetDistInAV();
        double distInWater = lightPath.GetDistInWater();
        fTransitTimes[i] = DU::Utility::Get()->GetGroupVelocity()
          .CalcByDistance(distInInnerAV, distInAV, distInWater,
                          RAT::util::WavelengthToEnergy(lambda * CLHEP::nm));
        fLEDIniDir[i] = lightPath.GetInitialLightVec();
      }
      // We have calculated all transittimes,
      // so no need to do this anymore for this run.
      fRunCalcTransitTimes = TRUE;
    }
    // Get the LED number, if we are in TELLIE case
    // Find out how many PMTs were hit in this event
    // Store this information to figure out occupancy later
    int numberCalPMTs = ev.GetCalPMTs().GetCount();
    fNHitCheck.push_back(numberCalPMTs);
    if (numberCalPMTs < fNHitsMin) {
      fNHitsMin = numberCalPMTs;
    }
    if (numberCalPMTs > fNHitsMax) {
      fNHitsMax = numberCalPMTs;
    }
    // Loop over all these hit PMTs and extract time and QHS QHL for each hit
    for (int involvPMT = 0; involvPMT < numberCalPMTs; involvPMT++) {
      const DS::PMTCal& pmtCal = ev.GetCalPMTs().GetPMT(involvPMT);
      const DS::PMTUncal& pmtUnCal = ev.GetUncalPMTs().GetPMT(involvPMT);
      int pmtID = pmtCal.GetID();
      // Set the position of this PMT
      pmtPos = pmtList.GetPosition(pmtID);
      // Do some checks before we decide that this is a hit we want to use:
      // 1) Check the ECA calibration status word
      if (pmtCal.GetStatus().GetULong64_t(0) != 0) {
        // Skip over this PMT hit since something was
        // fishy with the ECA calibration
        continue;
      }
      // 2) Is this a normal PMT?
      if (pmtList.GetType(pmtID) != 1) {
        // Skip over this PMT hit since it involves an abnormal PMT
        continue;
      }
      // 3) Viewing angle
      // Refracted angle
      double injectThetaTrue = TMath::RadToDeg()
        * TMath::ACos(fLEDIniDir[pmtID].Dot(ledDir)
                      / (fLEDIniDir[pmtID].Mag()
                         * ledDir.Mag()));
      // Fill the QHS and QHL vectors for this PMT
      fQHSCal->Fill(pmtCal.GetQHS());
      fQHLCal->Fill(pmtCal.GetQHL());
      fQHSUnCal->Fill(pmtUnCal.GetQHS());
      fQHLUnCal->Fill(pmtUnCal.GetQHL());
      fTimeCal->Fill(pmtCal.GetTime());
      fTimeUnCal->Fill(pmtUnCal.GetTime());
      fTimeResidual->Fill(pmtCal.GetTime() - fTransitTimes[pmtID]);
      // Count the hits per PMT
      fHitCount[pmtID]++;
      // Count the hits per PMT in the central region
      if (std::abs(injectThetaTrue) < fSelectAngle) {
        fHitCountCentral[pmtID]++;
      }
    }
    return OK;
  }


  int PCATellieMonitorProc::Save1DHisto(TH1F *hist,
                                        const std::string &axisTitle)
  {
    TImage *img = TImage::Create();
    TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);
    char imagename[30];
    snprintf(imagename, sizeof(imagename), "%s-%i.png",
             hist->GetTitle(), fRunID);
    c1->cd();
    c1->Clear();
    c1->SetLogy();
    hist->GetXaxis()->SetTitle(axisTitle.c_str());
    hist->GetYaxis()->SetTitle("Counts");
    hist->SetAxisRange(hist->GetBinCenter(hist->FindFirstBinAbove(0)) -5,
                       hist->GetBinCenter(hist->FindLastBinAbove(0)) +5, "X");
    hist->SetTitle("");
    hist->SetStats(0);
    hist->Draw();
    img->Clear();
    img->FromPad(c1);
    img->WriteImage(imagename);
    if (c1) c1->Close();
    img->Delete();
    return 1;
  }

  int PCATellieMonitorProc::SaveTGraph2D(TGraph2D *graph,
                                         const std::string &title)
  {
    TImage *img = TImage::Create();
    TView3D *view = new TView3D();
    TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);
    char imagename[30];
    snprintf(imagename, sizeof(imagename), "%s-%i.png", title.c_str(), fRunID);
    c1->cd();
    c1->Clear();
    c1->SetLogy(0);
    c1->SetLogz(0);
    graph->SetMarkerStyle(8);
    graph->SetMarkerSize(0.5);
    graph->GetXaxis()->SetLabelSize(0.03);
    graph->GetYaxis()->SetLabelSize(0.03);
    graph->GetXaxis()->SetRangeUser(0.0, 1.0);
    graph->GetYaxis()->SetRangeUser(0.0, 0.45);
    graph->SetTitle("");
    graph->Draw("ZCOLPCOL");
    view->RotateView(270, 0);
    img->Clear();
    img->FromPad(c1);
    img->WriteImage(imagename);
    if (c1) c1->Close();
    img->Delete();
    return 1;
  }
}   // namespace RAT