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

// ROOT stuff
#include <TCanvas.h>
#include <TFile.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TImage.h>
#include <TMath.h>
#include <TPaveStats.h>
#include <TStyle.h>
#include <TVector.h>
#include <TView3D.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/BinLog.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", 4096, 0, 4096);
    fQHLUnCal = new TH1F("QHLuncal", "QHLUnCal", 4096, 0, 4096);
    fTimeCal = new TH1F("Timecal", "TimeCal", 4000, -1000, 3000);
    fTimeUnCal = new TH1F("Timeuncal", "TimeUnCal", 4096, 0, 4096);
    fTimeResidual = new TH1F("Timeresidual", "TimeResidual", 4000, -1000, 3000);
    fTimeDiff = new TH1F("Timediff", "TimeDiff", 1000, 0, 1e9);
  }

  void PCATellieMonitorProc::BeginOfRun(DS::Run&)
  {

    // Changes by Martti, taken from DQTellieProc
    DBLinkPtr runInfo = DB::Get()->GetLink( "TELLIE_RUN","");
    json::Value subrunInfo = runInfo->GetJSON("sub_run_info");
    json::Value singleSubrunInfo = subrunInfo.getIndex(0);
    string fibre = singleSubrunInfo.getMember("fibre").getString();

    // 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 (fInputFibre) 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 fibre from the database
    /*
    std::ostringstream oss;
    oss.str("");
    // Take care of zero-padding
    oss << "FT" << std::setw(3) << std::setfill('0') << fInputFibre;
    // hard-coded to position 'A' at this time.
    // FIXME: in future, we need to find out which fibres are used
    // from an 'installed/used' entry in the database, not from user-input
    oss << "A";
    */
    try {
      //ledInfo =  DB::Get()->GetLink("FIBRE", oss.str());
      ledInfo =  DB::Get()->GetLink("FIBRE", fibre);
      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(): Fibre " << oss.str()
      info << "PCATellieMonitorProc(): Fibre " << fibre
           << " does not exist." << endl;
      // We should not proceed if we were unable to retrieve the fibreposition
      // 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)
        fInputFibre = value;
      else
        throw ParamInvalid(param, "value must be between 0 and 115");
    } else {
      throw ParamUnknown(param);
    }
  }

  void PCATellieMonitorProc::Smooth(TGraph* g, TGraph* h, int width)
  {
    int n = h->GetN();
    double* y = h->GetY();
    double sum;
    int cnt, min, max, b;
    for (int i=0; i<n; i++) {
      if (width==0) {
        g->SetPoint(i,i,y[i]);
        continue;
      }
      min = (i<width) ? 0 : i-width;
      max = (i>n-width) ? n : i+width;
      sum = 0;
      cnt = 0;
      for (b=min; b<max; b++) {
        sum += y[b];
        cnt++;
      }
      g->SetPoint(i,i,sum/cnt);
    }
  }

  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]");
    PCATellieMonitorProc::Save1DHisto(fTimeDiff, "Time_diff [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 min/max number of hits per PMT.
    int maxHits = -9000;
    int minHits = 9000;

    // More histograms!
    TH1F *hitsPerPMT = new TH1F("HitsPerPMT", "HitsPerPMT", 200, 0, fEventCounter);
    TH1F *hitsPerPMTAll = new TH1F("HitsPerPMTAll", "HitsPerPMTAll", 200, 0, fEventCounter);
    BinLog(hitsPerPMT->GetXaxis(), 0.5); // minimum for log scale
    BinLog(hitsPerPMTAll->GetXaxis(), 0.5); // minimum for log scale
    TH1F *occPerPMT = new TH1F("OccPerPMT", "OccPerPMT", 1000, 0, 1);

    // TGraph2D for the detector flatmaps
    TGraph2D *coverage = new TGraph2D();
    TGraph2D *coverageSelect = new TGraph2D();

    // Get PMT info 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);

    // Loop over all channels, fill the flatmaps and do some bookkeeping
    int graphIndex = 0;
    for (int i = 0; i < fTotPMTs; i++) {

      // Care only about normal (and HQE) PMTs for FlatMap
      if (pmtList.GetType(i)!=1 && pmtList.GetType(i)!=7) continue;

      // Use the FlatMap function to get the 2D position
      TVector3 pmtPos = pmtList.GetPosition(i);
      TVector2 pmtPosFlat;
      RAT::SphereToIcosahedron(pmtPos, pmtPosFlat, 1, 2.12);

      // Check for NAN or zero positions
      if (pmtPosFlat.X() != pmtPosFlat.X() || pmtPosFlat.X() == 0) {
        warn << "*** Warning *** PMT #" << i << " has position ( ";
        warn << pmtPos.X() << " | " << pmtPos.Y() << " | " << pmtPos.Z();
        warn << " ), can't be projected to FlatMap!" << newline;
        continue;
      }

      // Check if PMT was online
      if (chs.IsTubeOnline(i)) {
        coverage->SetPoint(graphIndex, pmtPosFlat.X(), pmtPosFlat.Y(), fHitCount[i]);
        coverageSelect->SetPoint(graphIndex, pmtPosFlat.X(), pmtPosFlat.Y(), fHitCountCentral[i]);
      } else { // PMT was offline, set negative value
        coverage->SetPoint(graphIndex, pmtPosFlat.X(), pmtPosFlat.Y(), -1);
        coverageSelect->SetPoint(graphIndex, pmtPosFlat.X(), pmtPosFlat.Y(), -1);
      }
      graphIndex++;

      // Fill histograms
      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", true);
    PCATellieMonitorProc::Save1DHisto(hitsPerPMTAll, "Number of hits", true);
    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]);
    }

    // Fit Gaussian to nhit distribution
    double max = nHits->GetMaximum();
    double lo = nHits->GetXaxis()->GetBinLowEdge(nHits->FindFirstBinAbove(0.2*max));
    double hi = nHits->GetXaxis()->GetBinUpEdge(nHits->FindLastBinAbove(0.2*max));
    TF1* gaussian = new TF1("gaussian", "gaus", lo, hi);
    nHits->Fit(gaussian, "Rq");

    // Fit Poisson distribution instead of Gaussian
    /*
    TF1* poisson = new TF1("poisson", "[0]*TMath::Poisson(x,[1])", max-50, max+50);
    poisson->SetParameter(0,5000);
    poisson->SetParameter(1,max);
    poisson->SetParName(0,"Normalisation");
    poisson->SetParName(1,"#mu");
    nHits->Fit(poisson, "Rq");
    */

    // Check nhit spread
    double avNhit = gaussian->GetParameter(1);
    double nhitSpread = gaussian->GetParameter(2);
    //double avNhit = poisson->GetParameter(1);
    //double nhitSpread = sqrt(poisson->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;
    }

    // Modify plot options and draw nhit distribution
    c1->cd();
    c1->Clear();
    c1->SetRightMargin(0.05);
    c1->SetTopMargin(0.05);
    gStyle->SetOptStat(0);
    gStyle->SetOptFit(1);
    nHits->SetAxisRange(0, nHits->FindLastBinAbove(0)+5);
    nHits->GetXaxis()->SetLabelSize(0.03);
    nHits->GetYaxis()->SetLabelSize(0.03);
    nHits->GetXaxis()->SetTitle("Nhit");
    nHits->GetYaxis()->SetTitle("Counts");
    nHits->GetXaxis()->SetTitleOffset(1.4);
    nHits->GetYaxis()->SetTitleOffset(1.4);
    nHits->SetTitle("");
    nHits->SetStats(1);
    nHits->Draw();
    gaussian->Draw("SAME");
    gPad->Modified(); gPad->Update(); // make sure it’s (re)drawn

    // Draw statistics box
    TPaveStats* st = (TPaveStats*)nHits->FindObject("stats");
    st->SetOptFit(1);
    gPad->Modified(); gPad->Update(); // make sure it’s (re)drawn
    st->SetX1NDC(0.70);
    st->SetX2NDC(0.95);
    st->SetY1NDC(0.73);
    st->SetY2NDC(0.95);

    // Save plot and close
    gPad->SetTicks(1,1);
    c1->Update();
    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];
    //TH1F* nhithist = new TH1F("nhithist","",numPoints,0,numPoints);
    for (size_t i = 0; i < numPoints; i++) {
      xarray[i] = i;
      yarray[i] = fNHitCheck[i];
      //nhithist->SetBinContent(i+1,fNHitCheck[i]);
    }
    TGraph *nhitGraph = new TGraph(numPoints, xarray, yarray);
    //TF1 *polynomial = new TF1("pol1", "pol1", 0, numPoints);
    //nhitGraph->Fit(polynomial, "Rq");
    // Smooth function
    TGraph* smooth = new TGraph();//("smooth","",numPoints,0,numPoints);
    Smooth(smooth,nhitGraph,1000);
    double meanS = smooth->GetMean(2);
    double *valS = smooth->GetY();
    double minS = 1e6;
    double maxS = -1e6;
    for (int i=0; i<smooth->GetN(); i++) {
      if (valS[i] > maxS) maxS = valS[i];
      if (valS[i] < minS) minS = valS[i];
    }
    double TOL = 0.2; // tolerance for fluctuations (20%)
    //if (TMath::Abs(polynomial->GetParameter(1)) - (polynomial->GetParError(1) / 3.0) < 0) {
    if (maxS/meanS > 1.+TOL || minS/meanS < 1.-TOL) {
      // 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);
    c1->SetRightMargin(0.05);
    c1->SetTopMargin(0.05);
    gStyle->SetOptStat(0);
    nhitGraph->GetXaxis()->SetTitle("Event number");
    nhitGraph->GetYaxis()->SetTitle("NHit");
    nhitGraph->GetXaxis()->SetTitleOffset(1.4);
    nhitGraph->GetYaxis()->SetTitleOffset(1.4);
    nhitGraph->SetTitle("");
    nhitGraph->Draw();
    smooth->SetLineWidth(2);
    smooth->SetLineColor(2);
    smooth->Draw("same");
    gPad->SetTicks(1,1);
    c1->Update();
    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", fInputFibre);
    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)
  {
    // Require EXTASYNC trigger - TODO also look for stolen triggers!
    int trig = ev.GetTrigType();
    if (trig != fTrigType) {
      return OK;
    }
    // TODO: fill fTimeDiff with trigger time differences here!
    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();
    const DU::PMTCalStatus& pmtStatus = DU::Utility::Get()->GetPMTCalStatus();
    // 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 operations 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 fibre 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();
      unsigned int status = pmtStatus.GetHitStatus(pmtCal);
      // 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 PMT status is good (as done in PMTCalSelector)
      if(status & (1<<pmtStatus.kCHSBit)) continue;
      if(status & (1<<pmtStatus.kECABit)) continue;
      //if(status & (1<<pmtStatus.kPCABit)) continue; // check for this after PCA
      if(status & (1<<pmtStatus.kXTalkBit)) continue;
      // 2) Is this a normal or HQE PMT?
      if (pmtList.GetType(pmtID)!=1 && pmtList.GetType(pmtID)!=7) continue;
      // 3) Viewing angle
      // Refracted angle
      double injectThetaTrue = TMath::RadToDeg()*fLEDIniDir[pmtID].Angle(ledDir);
      // 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,
                                        const bool logX)
  {
    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();
    if (logX) c1->SetLogx();
    c1->SetLogy();
    c1->SetRightMargin(0.05);
    c1->SetTopMargin(0.05);
    hist->GetXaxis()->SetTitle(axisTitle.c_str());
    hist->GetYaxis()->SetTitle("Counts");
    hist->GetXaxis()->SetTitleOffset(1.4);
    hist->GetYaxis()->SetTitleOffset(1.4);
    double minval = hist->GetBinCenter(hist->FindFirstBinAbove(0)-5);
    double maxval = hist->GetBinCenter(hist->FindLastBinAbove(0)+5);
    while (logX && maxval/minval<30.) { // expand range for log scale
      //std::cout << "Min = " << minval << ", Max = " << maxval << std::endl;
      if (minval <= 0.) { // protection for loop
        minval = maxval*1e-4;
        break;
      }
      minval /= 2;
      maxval *= 2;
    }
    hist->SetAxisRange(minval,maxval,"X");
    hist->SetTitle("");
    hist->SetStats(0);
    hist->Draw();
    gPad->SetTicks(1,1);
    c1->Update();
    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)
  {
    const int NCOL = 50;

    int NPTS = graph->GetN();
    double *x = graph->GetX();
    double *y = graph->GetY();
    double *z = graph->GetZ();

    double MAXVAL = graph->GetZmax();
    double MINVAL = MAXVAL / 1e3;

    std::vector<TGraph*> g(NCOL+1, NULL);
    for (int i=0; i<NPTS; i++) {
      int step;

      // Get correct "step" for colour scale
      if (z[i] < 0) step=0; // offline PMT
      else if (z[i] < MINVAL) step=1; // low rate (or bad DQ)
      // linear colour scale
      //else step = 1 + (int)TMath::Ceil(z[i]/(1.*MAXVAL/(NCOL-1)));
      // logarithmic colour scale
      else step = 1 + (int)TMath::Ceil( (log10(z[i])-log10(MINVAL)) /
                                        ((log10(MAXVAL)-log10(MINVAL))/(NCOL-1)) );

      // Add PMT to correct graph in vector
      if (!g[step]) g[step] = new TGraph();
      g[step]->SetPoint(g[step]->GetN(),x[i],y[i]);
    }

    // Set color scale to include offline PMTs
    int colors[NCOL+1];
    for (int i=0; i<=NCOL; i++) {
      if (i==0) colors[i] = 16; // grey
      else      colors[i] = (int)(50.+50./NCOL*i);
    }
    gStyle->SetNumberContours(NCOL);
    gStyle->SetPalette(NCOL+1,colors);

    // Create image and canvas
    TImage *img = TImage::Create();
    TCanvas *c1 = new TCanvas("c1", "c1", 1500, 600);
    char imagename[30];
    snprintf(imagename, sizeof(imagename), "%s-%i.png", title.c_str(), fRunID);
    c1->cd();
    c1->SetFrameLineColor(0);
    c1->SetFrameBorderSize(0);
    c1->SetFrameBorderMode(0);
    c1->SetLeftMargin(0.01);
    c1->SetRightMargin(0.09);
    c1->SetTopMargin(0.02);
    c1->SetBottomMargin(0.02);
    c1->SetLogz();

    // Draw frame with z-axis colour scale
    TH2F h;
    h.SetBins(10,0,1,10,0,0.45);
    h.SetStats(0);
    h.Fill(-999,-999);
    h.Draw("colza");
    h.GetZaxis()->SetRangeUser(MINVAL,MAXVAL);

    // Draw graphs
    for (int i=0; i<=NCOL; i++) {
      if (!g[i]) continue;
      g[i]->SetMarkerStyle(8);
      g[i]->SetMarkerSize(0.5);
      g[i]->SetMarkerColor(colors[i]);
      g[i]->Draw("P same");
    }

    // Save graphs to file
    c1->SetTicks(0,0);
    c1->Update();
    img->Clear();
    img->FromPad(c1);
    img->WriteImage(imagename);
    if (c1) c1->Close();
    img->Delete();
    return 1;
  }

}   // namespace RAT