// See DQTellieProc.hh for details
#include <RAT/Processor.hh>
#include <RAT/DataQualityProc.hh>
#include <RAT/DQTellieProc.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/EV.hh>
#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/DBLink.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/json.hh>
#include <RAT/BinLog.hh>

#include <sstream>

#include <TCanvas.h>
#include <TGraph2D.h>
#include <TF1.h>
#include <TF2.h>
#include <TH1I.h>
#include <TH2D.h>
#include <TLegend.h>
#include <TMarker.h>
#include <TMath.h>
#include <TPad.h>
#include <TSpectrum.h>
#include <TStyle.h>
#include <TVectorD.h>
#include <TVector2.h>
#include <TVector3.h>

using namespace std;
using namespace json;

namespace RAT {

  DQTellieProc::DQTellieProc() : DataQualityProc( "dqtellieproc")
  {
    fDQChecks = NULL;
    fTrigCheckThresh = 0;
    fPulseSeparationThresh = 0;
  }

  DQTellieProc::~DQTellieProc()
  {

  }

  void DQTellieProc::BeginOfRun( DS::Run& run )
  {
    detail << "Starting TELLIE DQ..." << newline;
    DataQualityProc::BeginOfRun( run );

    // Get variables from DB links
    fDQChecks = DB::Get()->GetLink( "DQCHECKS", "tellie" );
    DBLinkPtr runInfo = DB::Get()->GetLink( "TELLIE_RUN","");

    fTrigCheckThresh = fDQChecks->GetD( "trigger_check_thresh" );
    fCriteria["trigger_check_thresh"] = fTrigCheckThresh;
    fPulseSeparationMin = fDQChecks->GetD( "pulse_separation_min" );
    fCriteria["pulse_separation_min"] = fPulseSeparationMin;
    fPulseSeparationMax = fDQChecks->GetD( "pulse_separation_max" );
    fCriteria["pulse_separation_max"] = fPulseSeparationMax;
    fTellieChipDelay  = fDQChecks->GetD( "tellie_chip_delay" );
    fCriteria["tellie_chip_delay"] = fTellieChipDelay;
    fPulseSeparationThresh = fDQChecks->GetD( "pulse_separation_thresh" );
    fCriteria["pulse_separation_thresh"] = fPulseSeparationThresh;
    fMaxNHit = fDQChecks->GetI( "max_nhit");
    fCriteria["max_nhit"] = fMaxNHit;
    fMaxAvgNHit = fDQChecks->GetD( "max_avg_nhit");
    fCriteria["max_avg_nhit"] = fMaxAvgNHit;
    fTACPromptHeightRatio = fDQChecks->GetD( "tac_peak_amplitude_ratio");
    fCriteria["tac_peak_amplitude_ratio"] = fTACPromptHeightRatio;
    fTACPromptHeightThresh = fDQChecks->GetD( "tac_prompt_amplitude_thresh");
    fCriteria["tac_prompt_amplitude_thresh"] = fTACPromptHeightThresh;
    fTACReflectedHeightThresh = fDQChecks->GetD( "tac_reflected_amplitude_thresh");
    fCriteria["tac_reflected_amplitude_thresh"] = fTACReflectedHeightThresh;
    fTACPromptSigmaThresh = fDQChecks->GetD( "tac_prompt_sigma_thresh");
    fCriteria["tac_prompt_sigma_thresh"] = fTACPromptSigmaThresh;
    fTACReflectedSigmaThresh = fDQChecks->GetD( "tac_reflected_sigma_thresh");
    fCriteria["tac_reflected_sigma_thresh"] = fTACReflectedSigmaThresh;
    fTACPromptPeakLowerLimit = fDQChecks->GetD("tac_prompt_peak_low");
    fCriteria["tac_prompt_peak_low"] = fTACPromptPeakLowerLimit;
    fTACPromptPeakUpperLimit =  fDQChecks->GetD("tac_prompt_peak_high");
    fCriteria["tac_prompt_peak_high"] = fTACPromptPeakUpperLimit;
    fTACEarlyOffsetLower = fDQChecks->GetD("tac_pre_peak_offset_low");
    fCriteria["tac_pre_peak_offset_low"] = fTACEarlyOffsetLower;
    fTACEarlyOffsetUpper = fDQChecks->GetD("tac_pre_peak_offset_high");
    fCriteria["tac_pre_peak_offset_high"] = fTACEarlyOffsetUpper;
    fTACLateOffsetLower = fDQChecks->GetD("tac_late_peak_offset_low");
    fCriteria["tac_late_peak_offset_low"] = fTACLateOffsetLower;
    fTACLateOffsetUpper = fDQChecks->GetD("tac_late_peak_offset_high");
    fCriteria["tac_late_peak_offset_high"] = fTACLateOffsetUpper;
    fTriggerBit = fDQChecks->GetI("trigger_bit");
    fCriteria["trigger_bit_checked"] = fTriggerBit;
    fMaxStolenTriggers = fDQChecks->GetI("stolen_triggers_max");
    fCriteria["stolen_triggers_max"] = fMaxStolenTriggers;

    //Get JSON array from run information file
    json::Value subrunInfo = runInfo->GetJSON("sub_run_info");
    unsigned int subrunArraySize = subrunInfo.getArraySize();
    for(unsigned int i=0; i<subrunArraySize; i++) {
      json::Value singleSubrunInfo = subrunInfo.getIndex(i);
      try {
        fSubRunNumbers.push_back(singleSubrunInfo.getMember("sub_run_number").getInteger());
      } catch (std::runtime_error) { // sub_run_info sometimes missing?
        fSubRunNumbers.push_back(i);
      }
      fFibres.push_back(singleSubrunInfo.getMember("fibre").getString());
      fRunMode.push_back(singleSubrunInfo.getMember("run_mode").getString());
      fTellieEvents.push_back(singleSubrunInfo.getMember("number_of_shots").getInteger());
      try {
        fTelliePulseRate.push_back(1000.0/singleSubrunInfo.getMember("pulse_separation").getReal());
      } catch (std::runtime_error) { // reading JSON tables is rubbish
        fTelliePulseRate.push_back(1000.0/singleSubrunInfo.getMember("pulse_separation").getInteger());
      }
    }

    fLastSubRunIndex = -1;

    //Intialize the vectors
    fTellieCount = std::vector<int>(fSubRunNumbers.size(),0);
    fStolenTriggerCount = std::vector<int>(fSubRunNumbers.size(),0);
    fFailPulseSeparationCount = std::vector<int> (fSubRunNumbers.size(),0);
    fAvgNHit = std::vector<double> (fSubRunNumbers.size(),0.0);
    fStolenEvents = std::vector<DS::EV> (fMaxStolenTriggers,DS::EV());
    fRecoveredEvents = std::vector<DS::EV> (fMaxStolenTriggers,DS::EV());

    fTH1DLcnHits  = std::vector<TH1D*> (fSubRunNumbers.size(),NULL);
    fTH1DTacCount  = std::vector<TH1D*> (fSubRunNumbers.size(),NULL);
    fTH1DTacCountReflected  = std::vector<TH1D*> (fSubRunNumbers.size(),NULL);
    fTH1DTacCountDirect  = std::vector<TH1D*> (fSubRunNumbers.size(),NULL);
    fTH1DCalibTime = std::vector<TH1D*> (fSubRunNumbers.size(),NULL);
    fTH1DPulseSeparation = std::vector<TH1D*> (fSubRunNumbers.size(),NULL);
    fTH1DPulseSeparationStolen = std::vector<TH1D*> (fSubRunNumbers.size(),NULL);
    fTH2DPsupProjection  = std::vector<TH2D*> (fSubRunNumbers.size(),NULL);
    fTH2DPsupOffline  = std::vector<TH2D*> (fSubRunNumbers.size(),NULL);
    fTH1INHits  = std::vector<TH1I*> (fSubRunNumbers.size(),NULL);

    fFibreFiringGuess = std::vector<std::string> (fSubRunNumbers.size(),"");
    fFibrePositions = std::vector<TVector3*> (fSubRunNumbers.size(), NULL);

    fPromptPeakADCCount = std::vector<double> (fSubRunNumbers.size(),0.0);
    fPrePeakADCCount = std::vector<double> (fSubRunNumbers.size(),0.0);
    fLatePeakADCCount = std::vector<double> (fSubRunNumbers.size(),0.0);

    fSubRunLengths = std::vector<double> (fSubRunNumbers.size(),0.0);
    fPulseSeparationEfficiency = std::vector<double> (fSubRunNumbers.size(),0.0);

    fPromptPeakADCCountCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fPromptPeakAmplitudeCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fPeakTimeSpacingCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fPeakNumberCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fTriggerCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fSubRunLengthsCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fPulseSeparationEfficiencyCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fAvgNHitCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fMaxNHitCheck = std::vector<int> (fSubRunNumbers.size(),0);
    fFailMaxNHitCount = std::vector<int> (fSubRunNumbers.size(),0);
    fCorrectFibreCheck = std::vector<int> (fSubRunNumbers.size(),0);

    fFirstEventTime = std::vector<DS::UniversalTime> (fSubRunNumbers.size(),DS::UniversalTime());
    fLastEventTime = std::vector<DS::UniversalTime> (fSubRunNumbers.size(),DS::UniversalTime());
    fPrevEXTATime = std::vector<DS::UniversalTime> (fSubRunNumbers.size(),DS::UniversalTime());
    fCurEventTime = std::vector<DS::UniversalTime> (fSubRunNumbers.size(),DS::UniversalTime());

    fFirstEventTimeSet = std::vector<bool> (fSubRunNumbers.size(),false);
    fFirstTellieEventTimeSet = std::vector<bool> (fSubRunNumbers.size(),false);

    fHitCount.resize(fSubRunNumbers.size(),std::vector<int>(0,0));

    fNPeaks = std::vector<int> (fSubRunNumbers.size(),3); //Pre, prompt and late pulse peaks

    //Initialise hit count
    for(unsigned int i=0; i<fSubRunNumbers.size(); i++) {
      // Resize for the given run
      fHitCount[i].resize(DU::Utility::Get()->GetPMTInfo().GetCount(), 0);
    }
    for(unsigned int i=0; i<fSubRunNumbers.size(); i++) {
      //Book Histograms
      stringstream stringRunID;
      stringstream namestream;
      string titleString;
      namestream << "_run_"<<fRunID;
      namestream << "_subrun_"<<fSubRunNumbers[i];
      stringRunID << "SNO+ Data Quality: TELLIE - run ";
      stringRunID << fRunID;
      stringRunID << " - subrun ";
      stringRunID << fSubRunNumbers[i];
      string baseTitle = stringRunID.str();
      titleString = "lcn_hits"+namestream.str();
      //For finding max hit PMT
      string title = baseTitle + " LCN Hits (Hot PMTs removed)";
      fTH1DLcnHits[i] = new TH1D(titleString.c_str(), title.c_str(),
          10000, 0, 10000);
      fTH1DLcnHits[i]->SetDirectory(0);
      //Timing checks
      title = baseTitle + " TAC Count";
      titleString = "tac_count"+namestream.str();
      fTH1DTacCount[i] = new TH1D(titleString.c_str(), title.c_str(),
          256,0,4096);
      fTH1DTacCount[i]->SetDirectory(0);
      //Direct TAC
      title = baseTitle + " TAC Count Direct";
      titleString = "tac_count_direct"+namestream.str();
      fTH1DTacCountDirect[i] = new TH1D(titleString.c_str(), title.c_str(),
          256,0,4096);
      fTH1DTacCountDirect[i]->SetDirectory(0);
      //Reflected TAC
      title = baseTitle + " TAC Count Reflected";
      titleString = "tac_count_reflected"+namestream.str();
      fTH1DTacCountReflected[i] = new TH1D(titleString.c_str(), title.c_str(),
          256,0,4096);
      fTH1DTacCountReflected[i]->SetDirectory(0);
      //For post run monitoring
      title = baseTitle + " Calibrated Time";
      titleString = "calibrated_time"+namestream.str();
      fTH1DCalibTime[i] = new TH1D(titleString.c_str(), title.c_str(),
          500, -500, 500);
      fTH1DCalibTime[i]->SetDirectory(0);
      // Trigger separation between EXTA events
      title = baseTitle + " Pulse separations";
      titleString = "pulse_separation"+namestream.str();
      double maxTrigSep = 0, binsPerMS = 0;
      if (fRunMode[i] == "Master") {
        maxTrigSep = fMaxStolenTriggers*(1000./fTelliePulseRate[i]+fTellieChipDelay);
        binsPerMS = 5e2; // 2 us bin width
      } else {
        maxTrigSep = fMaxStolenTriggers*(1000./fTelliePulseRate[i]);
        binsPerMS = 2e4; // 50 ns bin width (to resolve 10 MHz clock)
      }
      fTH1DPulseSeparation[i] = new TH1D(titleString.c_str(), title.c_str(),
          binsPerMS*maxTrigSep, -0.5/binsPerMS, maxTrigSep-0.5/binsPerMS);
      fTH1DPulseSeparation[i]->SetDirectory(0);
      // Trigger separation between other events suspected to have stolen EXTA trigger
      title = baseTitle + " Stolen pulse separations";
      titleString = "stolen_pulse_separation"+namestream.str();
      fTH1DPulseSeparationStolen[i] = new TH1D(titleString.c_str(), title.c_str(),
          binsPerMS*maxTrigSep, -0.5/binsPerMS, maxTrigSep-0.5/binsPerMS);
      fTH1DPulseSeparationStolen[i]->SetDirectory(0);
      //For post run monitoring
      title = baseTitle + " PSUP Hit Map";
      titleString = "psup_projection"+namestream.str();
      fTH2DPsupProjection[i] = new TH2D(titleString.c_str(), title.c_str(),
          400, 0, 1, 200, 0, 1);
      fTH2DPsupProjection[i]->SetDirectory(0);
      fTH2DPsupOffline[i] = new TH2D(titleString.c_str(), "",
          400, 0, 1, 200, 0, 1);
      fTH2DPsupOffline[i]->SetDirectory(0);
      title = baseTitle + " NHits ";
      titleString = "nhits"+namestream.str();
      fTH1INHits[i] = new TH1I(titleString.c_str(), title.c_str(),
          1000, 0, 1000);
      fTH1INHits[i]->SetDirectory(0);
      stringstream fibNum;
      if (fFibres[i].size() == 6)
      {
        fibNum << fFibres[i];
        debug << "fibNum" << fibNum.str() << newline;
      }
      else if(*fFibres[i].rbegin() == 'A')
      {
        fibNum << "FT" << setw(4) << setfill('0') << fFibres[i];
        fFibres[i] = fibNum.str();
        debug  << "fibNum2" << fibNum.str() << newline;
      }
      else
      {
        fibNum << "FT" << setw(3) << setfill('0') << fFibres[i] << "A";
        fFibres[i] = fibNum.str();
        debug  << "fibNum3" << fibNum.str() << newline;
      }
      DBLinkPtr fibre = DB::Get()->GetLink("FIBRE", fibNum.str());
      fFibrePositions[i] = new TVector3(0., 0., 0.);
      fFibrePositions[i]->SetX(fibre->GetD("x"));
      fFibrePositions[i]->SetY(fibre->GetD("y"));
      fFibrePositions[i]->SetZ(fibre->GetD("z"));
    }
  }

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

  Processor::Result DQTellieProc::Event( DS::Entry& ent, DS::EV& ev )
  {
    // Get the subrun number and find its index within vector of subruns
    int subRunID = ent.GetSubRunID();
    int subRunIndex = -1;
    for(unsigned int i=0; i<fSubRunNumbers.size(); i++) {
      if (subRunID != fSubRunNumbers[i]) continue;
      subRunIndex = i;
      break;
    }

    // Subrun number not found in couchDB subruns
    if (subRunIndex == -1) {
      detail << "Event with subrun " << subRunID;
      detail << " is not in CouchDB subrun range: [";
      detail << fSubRunNumbers.front() << ",";
      detail << fSubRunNumbers.back() << "]" << newline;
      detail << "Skipping this event." << newline;
      return OK;
    }

    const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
    // Check 1 - Check trigger mask
    bool passedTriggerCheck = false;
    if (ev.GetTrigType() & (1<<fTriggerBit) ) {
      // EXTA trigger as expected
      passedTriggerCheck = true;
      fTellieCount[subRunIndex]++;
    }

    // Check 2 - Check run time make sure we aren't getting times from orphans
    if( fFirstEventTimeSet[subRunIndex] == false && ev.GetGTID() !=0) {
      fFirstEventTime[subRunIndex] = ev.GetUniversalTime();
      fFirstEventTimeSet[subRunIndex] = true;
    } else if(ev.GetGTID() != 0) {
      fLastEventTime[subRunIndex] = ev.GetUniversalTime();
    }

    //Check 3 - Check pulse separation
    bool passedSeparationCheck = false;
    bool multipleSeparationCheck = false;
    if( fFirstTellieEventTimeSet[subRunIndex] == false && ev.GetGTID() != 0 ) {
      fPrevEXTATime[subRunIndex] = ev.GetUniversalTime();
      fFirstTellieEventTimeSet[subRunIndex] = true;
    } else if(ev.GetGTID() != 0) {

      // trigger separation = nominal + allowed offset + chip delay (Master mode)
      double nominalPulseSeparation = 1.0/fTelliePulseRate[subRunIndex];
      double lowerPulseSeparation = nominalPulseSeparation + fPulseSeparationMin;
      double upperPulseSeparation = nominalPulseSeparation + fPulseSeparationMax;
      if (fRunMode[subRunIndex] == "Master") { // add internal PIC chip delay
        // allowed trigger separation range must be much less strict in Master mode
        lowerPulseSeparation += fTellieChipDelay + 99.*fPulseSeparationMin;
        upperPulseSeparation += fTellieChipDelay + 99.*fPulseSeparationMax;
        nominalPulseSeparation = (lowerPulseSeparation+upperPulseSeparation)/2.;
      }

      // calculate separation between this event and previous EXTA event
      fCurEventTime[subRunIndex] = fLastEventTime[subRunIndex];
      DS::UniversalTime utSeparation  = fCurEventTime[subRunIndex] - fPrevEXTATime[subRunIndex];
      double pulseSeparation = utSeparation.GetSeconds()+(utSeparation.GetNanoSeconds()*1e-9);
      passedSeparationCheck = (pulseSeparation > lowerPulseSeparation && pulseSeparation < upperPulseSeparation);
      multipleSeparationCheck = (pulseSeparation > 2*lowerPulseSeparation && pulseSeparation < fMaxStolenTriggers*upperPulseSeparation);

      if ( passedTriggerCheck ) { // EXTA trigger
        if ( !passedSeparationCheck) { // EXTA trigger separation not as expected
          if (subRunIndex != fLastSubRunIndex) { // first event per subrun
            fLastSubRunIndex = subRunIndex; // update previous subrun index
          } else if (multipleSeparationCheck) { // skipped multiple events?
            int nSkipped = (int)round(pulseSeparation/nominalPulseSeparation)-1;
            int nRecovered = 0;
            info << "INFO - Event " << Form("%8d", ev.GetGTID());
            info << " has pulse separation " << Form("%6.3f ms", 1000.*pulseSeparation);
            info << ", looking for " << nSkipped;
            info << " potentially stolen events..." << newline;
            if (fStolenEvents.size() == 0) {
              info << " ==> none found!" << newline;
              warn << "WARNING - Event " << Form("%8d", ev.GetGTID());
              warn << " has bad pulse separation " << Form("%6.3f ms", 1000.*pulseSeparation);
              warn << " (min/max " << 1000.*lowerPulseSeparation << "/";
              warn << 1000.*upperPulseSeparation << " ms)";
              warn << ", trigger word:";
              for (int t=0; t<=26; t++) {
                if (ev.GetTrigType() & (1<<t)) {
                  warn << " " << DU::TrigBits::GetBitName(t);
                }
              }
              warn << newline;
            } else {

              // Loop over potentially stolen events and tag those with highest NHit
              // FIXME - should check for most nhits within fibre aperture instead
              std::vector <int> stolenGTIDs;
              stolenGTIDs.clear();
              while ((int)stolenGTIDs.size() < nSkipped) {
                if(stolenGTIDs.size() == fStolenEvents.size()) break;
                int highNhit = -999;
                int highGTID = -999;
                for (std::vector<DS::EV>::iterator it = fStolenEvents.begin(); it != fStolenEvents.end(); ++it) {
                  DS::EV& thisEV = *it;
                  // skip event if it was already tagged as stolen trigger
                  if (std::find(stolenGTIDs.begin(), stolenGTIDs.end(), thisEV.GetGTID()) != stolenGTIDs.end()) continue;
                  if ((int)thisEV.GetNhits() > highNhit) {
                    highGTID = thisEV.GetGTID();
                    highNhit = thisEV.GetNhits();
                  }
                }
                stolenGTIDs.push_back(highGTID); // tag GTID as stolen trigger
              }

              // Recover tagged events
              fRecoveredEvents.clear();
              for (std::vector<DS::EV>::iterator it = fStolenEvents.begin(); it != fStolenEvents.end(); ++it) {
                DS::EV& thisEV = *it;
                // Skip event if it was not tagged as stolen trigger
                if (std::find(stolenGTIDs.begin(), stolenGTIDs.end(), thisEV.GetGTID()) == stolenGTIDs.end()) continue;
                // Skip event if the timing is too close to the current EXTA pulse
                DS::UniversalTime itSeparation = thisEV.GetUniversalTime() - fPrevEXTATime[subRunIndex];
                double thisSeparation = itSeparation.GetSeconds()+(itSeparation.GetNanoSeconds()*1e-9);
                if((int)round(thisSeparation/nominalPulseSeparation)-1 >= nSkipped) continue;

                // FIXME - these events should be tagged/recovered in separate event loops
                info << " ==> Event " << Form("%8d", thisEV.GetGTID());
                info << " has pulse separation " << Form("%6.3f ms", 1000.*thisSeparation);
                info << ", trigger word:";
                for (int t=0; t<=26; t++) {
                  if (thisEV.GetTrigType() & (1<<t)) {
                    info << " " << DU::TrigBits::GetBitName(t);
                  }
                }
                info << ", " << thisEV.GetNhits() << " nhits";
                info << " - stolen TELLIE event! Recovering..." << newline;

                // Recover stolen event and fill histograms
                fRecoveredEvents.push_back(thisEV); // recover event
                fTH1DPulseSeparation[subRunIndex]->Fill(1000.*thisSeparation); // ms
                fTH1DPulseSeparationStolen[subRunIndex]->Fill(1000.*thisSeparation); // ms
                fStolenTriggerCount[subRunIndex]++;

                // Update previous event time and pulse separation
                fPrevEXTATime[subRunIndex] = thisEV.GetUniversalTime(); // update previous EXTA time
                utSeparation  = fCurEventTime[subRunIndex] - fPrevEXTATime[subRunIndex];
                pulseSeparation = utSeparation.GetSeconds()+(utSeparation.GetNanoSeconds()*1e-9);
              }

            }
            nRecovered = fRecoveredEvents.size();
            info << "...successfully recovered " << nRecovered << " events" << newline;
            fFailPulseSeparationCount[subRunIndex] += nSkipped-nRecovered;
          } else { // EXTA event out of sync
            warn << "WARNING - Event " << Form("%8d", ev.GetGTID());
            warn << " has bad pulse separation " << Form("%6.3f ms",1000.*pulseSeparation);
            warn << " (min/max " << 1000.*lowerPulseSeparation << "/";
            warn << 1000.*upperPulseSeparation << " ms)";
            warn << ", trigger word:";
            for (int t=0; t<=26; t++) {
              if (ev.GetTrigType() & (1<<t)) {
                warn << " " << DU::TrigBits::GetBitName(t);
              }
            }
            warn << newline;
            fFailPulseSeparationCount[subRunIndex]++;
          }
        }
        // Delete potentially stolen events and fill histogram
        fStolenEvents.clear();
        fTH1DPulseSeparation[subRunIndex]->Fill(1000.*pulseSeparation); // ms
        fPrevEXTATime[subRunIndex] = fCurEventTime[subRunIndex]; // update previous EXTA time
      } else { // not EXTA trigger
        if ( passedSeparationCheck ) { // nominal EXTA trigger separation
          // potential TELLIE event stolen by other triggers (e.g. ESumH)
          fStolenEvents.push_back(ev);
        } else if (multipleSeparationCheck) { // skipped multiple events?
          for(int m=2; m<=fMaxStolenTriggers; m++) {
            if (pulseSeparation > m*lowerPulseSeparation && pulseSeparation < m*upperPulseSeparation) {
              fStolenEvents.push_back(ev);
            }
          }
        } else { // trigger separation not as expected
          // normal non-TELLIE event, ignore!
          return OK;
        }
      }
    }

    // Continue with EXTA events only (FIXME - ignores recovered events for now)
    if ( !passedTriggerCheck ) return OK;

    //Checks 4-5 - NHit Checks
    const int NHit = ev.GetNhits();
    if (NHit > fMaxNHit) fFailMaxNHitCount[subRunIndex]++;
    fTH1INHits[subRunIndex]->Fill(NHit);
    fAvgNHit[subRunIndex] += static_cast<double>(NHit);

    //Check 6 - Fibre Check
    for(unsigned int ipmt=0;ipmt<ev.GetCalPMTs().GetCount();ipmt++) {
      RAT::DS::PMTCal& pmt = ev.GetCalPMTs().GetPMT(ipmt);
      int lcn = pmt.GetID();
      fHitCount[subRunIndex][lcn] = fHitCount[subRunIndex][lcn]+1;
      fTH1DLcnHits[subRunIndex]->Fill(lcn);
      fTH1DCalibTime[subRunIndex]->Fill(pmt.GetTime());
      TVector3 pmtPos = pmtInfo.GetPosition(lcn);
      if (pmtPos.Mag() == 0) {
        warn << "WARNING: DQTellieProc::Event pmt lcn = " << lcn;
        warn << " is at 0,0,0! Skipping" << newline;
        continue;
      }
      int face;
      TVector2 projection = IcosProject(pmtPos,face);
      fTH2DPsupProjection[subRunIndex]->Fill(projection.X(),projection.Y());
    }

    //Checks 7-10 - Timing Checks
    for(unsigned int ipmt=0;ipmt<ev.GetUncalPMTs().GetCount();ipmt++) {
      RAT::DS::PMTUncal& uncalibratedPmt = ev.GetUncalPMTs().GetPMT(ipmt);
      fTH1DTacCount[subRunIndex]->Fill(uncalibratedPmt.GetTime());
      TVector3 pmtPosUnit = pmtInfo.GetPosition(uncalibratedPmt.GetID()).Unit();
      double unitDot = pmtPosUnit.Dot(fFibrePositions[subRunIndex]->Unit());
      if (unitDot < TMath::Cos(165.*TMath::Pi()/180.)) {
        fTH1DTacCountDirect[subRunIndex]->Fill(uncalibratedPmt.GetTime());
      } else if (unitDot > TMath::Cos(15.*TMath::Pi()/180.)) {
        fTH1DTacCountReflected[subRunIndex]->Fill(uncalibratedPmt.GetTime());
      }
    }

    return OK;
  }

  void DQTellieProc::PlotFibreCheck(const TVector2& directProject,const TVector2& reflectedProject, const std::vector<TVector2>& directCircle, const std::vector<TVector2>& reflectedCircle, const TVector2& likelyFibreProject, const TVector2& actualFibreProject, const string& assumedFibre, int subRunIndex) {

    // Create canvas and pads
    TCanvas *C1 = new TCanvas("C1","",1000,600);
    TPad * histoPad = new TPad("histoPad","histoPad",0.0,0.2,1.0,1.0);
    TPad * legendPad = new TPad("legendPad","legendPad",0.0,0.0,1.0,0.2);
    fTH2DPsupProjection[subRunIndex]->SetStats(0);
    C1->cd();
    histoPad->cd();
    histoPad->SetFrameLineColor(0);
    histoPad->SetLeftMargin(0.01);
    histoPad->SetBottomMargin(0.01);
    histoPad->SetLogz();

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

    // Plot flat map projection of PSUP, including offline PMTs
    float MINOCC = 1.e-4;
    float MAXOCC = 1.e-1;
    fTH2DPsupProjection[subRunIndex]->Scale(1./fTellieEvents[subRunIndex]);
    fTH2DPsupProjection[subRunIndex]->Draw("colza");
    fTH2DPsupProjection[subRunIndex]->GetZaxis()->SetRangeUser(MINOCC,MAXOCC);
    fTH2DPsupOffline[subRunIndex]->Scale(MINOCC);
    fTH2DPsupOffline[subRunIndex]->Draw("col same");
    fTH2DPsupOffline[subRunIndex]->GetZaxis()->SetRangeUser(MINOCC,MAXOCC);

    // Draw markers and contours for fit positions
    TMarker directSpot;
    TMarker reflectedSpot;
    TMarker likelyFibre;
    TMarker actualFibre;
    directSpot.SetMarkerStyle(2);
    reflectedSpot.SetMarkerStyle(3);
    likelyFibre.SetMarkerStyle(7);
    actualFibre.SetMarkerStyle(4);
    TGraph directContour, reflectedContour;
    for (unsigned int i=0; i<directCircle.size(); i++)
      directContour.SetPoint(i,directCircle[i].X(),directCircle[i].Y());
    for (unsigned int i=0; i<reflectedCircle.size(); i++)
      reflectedContour.SetPoint(i,reflectedCircle[i].X(),reflectedCircle[i].Y());
    directSpot.DrawMarker(directProject.X(),directProject.Y());
    reflectedSpot.DrawMarker(reflectedProject.X(),reflectedProject.Y());
    likelyFibre.DrawMarker(likelyFibreProject.X(),likelyFibreProject.Y());
    actualFibre.DrawMarker(actualFibreProject.X(),actualFibreProject.Y());
    directContour.Draw("P same");
    reflectedContour.Draw("P same");

    // Draw legend in separate pad
    TLegend * leg = new TLegend(0.1,0.05,0.9,1.0);
    string assumed = Form("Assumed fibre position (%s) ", assumedFibre.c_str());
    string firing = Form("Firing fibre position (%s) ", fFibres[subRunIndex].c_str());
    leg->AddEntry(&directSpot,"Direct light fit position","p");
    leg->AddEntry(&reflectedSpot,"Reflected light fit position","p");
    leg->AddEntry(&likelyFibre,assumed.c_str(),"p");
    leg->AddEntry(&actualFibre,firing.c_str(),"p");
    legendPad->cd();
    leg->SetBorderSize(0);
    leg->Draw();

    // Draw pads and save graphs to file
    C1->cd();
    legendPad->Draw();
    histoPad->Draw();
    histoPad->SetTicks(0,0);
    C1->Update();
    string filename = fTH2DPsupProjection[subRunIndex]->GetName();
    filename = filename+".png";
    stringstream namestream;
    namestream << "PSUP Hit Map subrun " << subRunIndex;
    C1->SetName(namestream.str().c_str());
    WriteToRoot(C1);
    C1->Print(filename.c_str());

    // Delete pointers
    delete leg;
    leg = NULL;
    delete histoPad;
    histoPad = NULL;
    delete legendPad;
    legendPad = NULL;
    delete C1;
    C1 = NULL;
  }

  // Check whether the correct fibre is firing
  void DQTellieProc::FibreCheck(int subRunIndex) {

    // Get PMT and fibre info
    const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
    RAT::DB *db = RAT::DB::Get();
    TVector3 fibrePos, fibreDir, lightPos;
    try {
      RAT::DBLinkPtr entry = db->GetLink("FIBRE",fFibres[subRunIndex]);
      fibrePos.SetXYZ(entry->GetD("x"), entry->GetD("y"), entry->GetD("z"));
      fibreDir.SetXYZ(entry->GetD("u"), entry->GetD("v"), entry->GetD("w"));
      lightPos = fibrePos + 2*fibrePos.Mag()*fibreDir; // projected light spot
    } catch ( DBNotFoundError &e ) {
      warn << "WARNING: DQTellieProc::FibreCheck " << fFibres[subRunIndex] << " does not exist. Skipping." << newline;
      return;
    }

    // Convert hit count to occupancy
    std::vector<float> occupancy;
    unsigned int NPMTS = fHitCount[subRunIndex].size();
    for(unsigned int i=0; i<NPMTS; i++) {
      occupancy.push_back((float)fHitCount[subRunIndex][i] / fTellieEvents[subRunIndex]);
    }

    // ********************************************************************
    // First iteration: estimate light spots based on PSUP face intensity
    // ********************************************************************

    // Get relative "heat" of each PSUP face (side of icosahedron containing PMT)
    TVector2 icospos;
    TVector3 pmtPos;
    float radiusPSUP=0, totalHeat=0;
    int nfacepmts[20] = {0}, nfacegoodpmts[20] = {0};
    TVector3 facecenter[20], faceweight[20];
    for(unsigned int id=0; id<NPMTS; id++) {
      if(pmtInfo.GetType(id) != DU::PMTInfo::NORMAL) continue;
      pmtPos = pmtInfo.GetPosition(id);
      if (pmtPos.Mag()==0) continue;
      radiusPSUP += pmtPos.Mag()*occupancy[id];
      totalHeat += occupancy[id];

      // Find out which PSUP face PMT is on
      int face;
      icospos = IcosProject(pmtPos,face);
      facecenter[face-1] += pmtPos;
      nfacepmts[face-1]++;

      // Calculate PSUP face heat
      faceweight[face-1] += pmtPos*occupancy[id];
      nfacegoodpmts[face-1]++;
    }
    radiusPSUP /= totalHeat;

    // Find "hottest" face in this subrun
    float faceheat[20]={0}, maxfaceheat=-1;
    for(unsigned int fc=0; fc<20; fc++) {
      if (nfacepmts[fc]==0) continue;
      if (nfacegoodpmts[fc]==0) continue;
      facecenter[fc] *= 1./nfacepmts[fc];   // all PMTs weighted equally
      faceweight[fc] *= 1./nfacegoodpmts[fc]; // good PMTS weighted by intensity
      faceheat[fc] = faceweight[fc].Mag();
      if (faceheat[fc]>maxfaceheat) {
        maxfaceheat=faceheat[fc];
      }
    }

    // Obtain a first guesstimate for light spot
    TVector3 bestguess(0,0,0);
    int hotfaces=0;
    for(unsigned int fc=0; fc<20; fc++) {
      // face is "hot" if intensity is >20% of hottest face
      if(faceheat[fc]<maxfaceheat/5.) continue;
      bestguess += faceweight[fc];
      hotfaces++;
    }
    if (hotfaces>0) {
      bestguess *= 1./hotfaces;
    }
    else {
      cerr<<"Something went wrong! No hot faces found."<<endl;
      fCorrectFibreCheck[subRunIndex] = 0;
      return;
    }
    bestguess.SetMag(radiusPSUP);
    TVector3 guess_dir =  bestguess;
    TVector3 guess_ref = -bestguess;
    // FIXME - improve guesstimate for reflected light spot?


    // ********************************************************************
    // Second iteration: take weighted average around estimated light spots
    // ********************************************************************

    // Max. opening angles for direct/reflected light (using fibre aperture)
    const float DIR_CONE = 48.;
    const float REF_CONE = 20.;
    TVector3 direct(0.,0.,0.);
    TVector3 reflected(0.,0.,0.);

    // Calculate weighted average for PMTs within cone of first estimate
    for(unsigned int id=0; id<NPMTS; id++) {

      // Check PMTs are ok
      if(pmtInfo.GetType(id) != DU::PMTInfo::NORMAL) continue;
      pmtPos = pmtInfo.GetPosition(id);
      if (pmtPos.Mag()==0) continue;

      // Check offline PMTs
      if (occupancy[id] == 0.) {
        //info << "PMT #" << id << " is offline or has zero occupancy!" << endl;
        int face;
        TVector2 offline = IcosProject(pmtPos,face);
        fTH2DPsupOffline[subRunIndex]->Fill(offline.X(),offline.Y());
      }

      // Reject bad occupancy PMTs for weight (FIXME - cuts should not be hardcoded)
      if (occupancy[id] < 0.001) continue;
      if (occupancy[id] > 0.1) continue;

      // Obtain weighted average for PMT positions within given cone
      if (pmtPos.Angle(guess_dir) < TMath::Pi()*DIR_CONE/180.) {
        direct += pmtPos*occupancy[id];
      } else if (pmtPos.Angle(guess_ref) < TMath::Pi()*REF_CONE/180.) {
        reflected += pmtPos*occupancy[id];
      }
    }

    // Safety check in case no PMTs were within cone
    if (direct.Mag() == 0) {
      warn << "*** WARNING *** No good PMTs in direct light cone! Setting to z-axis." << endl;
      direct.SetXYZ(0,0,1);
    }
    if (reflected.Mag() == 0) {
      warn << "*** WARNING *** No good PMTs in reflected light cone! Setting to z-axis." << endl;
      reflected.SetXYZ(0,0,1);
    }

    // Scale estimated light spot to average PSUP radius
    direct.SetMag(radiusPSUP);
    reflected.SetMag(radiusPSUP);

    // Fill graphs with new view, centered on weighted light spot
    TGraph2D *gDir = new TGraph2D();  // 2D graph (for 3D view and fit)
    TGraph2D *gRef = new TGraph2D();
    FillHemisphere(direct, occupancy, gDir);
    FillHemisphere(reflected, occupancy, gRef);

    // Get rotation angles for rotating view over weighted light spots
    double rot_Z1, rot_X1, rot_Z2, rot_X2;
    GetRotationAngles(direct, rot_Z1, rot_X1);
    GetRotationAngles(reflected, rot_Z2, rot_X2);


    // ********************************************************************
    // Third iteration: perform 2D Gaussian fit around weighted light spots
    // ********************************************************************

    // Fit 2D graphs with Gaussian surface
    const int NPARS = 4;
    double parDir[NPARS], parRef[NPARS];
    debug << "Performing 2D Gaussian fits..." << endl;
    FitLightSpot(gDir,radiusPSUP/1e3,DIR_CONE,parDir); // units: dist [m], ang [deg]
    FitLightSpot(gRef,radiusPSUP/1e3,REF_CONE,parRef);

    // Output fit results (verbose mode)
    debug << "FIT RESULTS:" << endl;
    debug << "- Direct light";
    for (int p=0; p<NPARS; p++) debug << " " << parDir[p];
    debug << endl;
    debug << "- Reflected light";
    for (int p=0; p<NPARS; p++) debug << " " << parRef[p];
    debug << endl;

    // Standard deviation (Gaussian sigma) translated to angle [deg]
    double sigmaAngDir = asin(1e3*parDir[3]/radiusPSUP)*180./TMath::Pi();
    double sigmaAngRef = asin(1e3*parRef[3]/radiusPSUP)*180./TMath::Pi();

    // Translate to point on PSUP sphere
    TVector3 directFit(1e3*parDir[1],1e3*parDir[2],sqrt(radiusPSUP*radiusPSUP-1e6*parDir[1]*parDir[1]-1e6*parDir[2]*parDir[2]));
    directFit.RotateY(rot_Z1);
    directFit.RotateZ(rot_X1);
    TVector3 reflectedFit(1e3*parRef[1],1e3*parRef[2],sqrt(radiusPSUP*radiusPSUP-1e6*parRef[1]*parRef[1]-1e6*parRef[2]*parRef[2]));
    reflectedFit.RotateY(rot_Z2);
    reflectedFit.RotateZ(rot_X2);

    // Save to new vectors
    TVector3 directLight = directFit;
    TVector3 reflectedLight = reflectedFit;

    // If fibre is affected by belly plates, revert to weighted method
    if (CheckBellyPlate(subRunIndex)) {
      info << "Direct light is affected by belly plates. ";
      info << "Using weighted average instead of Gaussian fit." << endl;
      directLight = direct;
      reflectedLight = reflected;
    }

    // Use previously used code to save parameters as usual
    TVector3 lightOrigin = radiusPSUP*(reflectedLight.Unit());
    stringstream namestream;
    namestream << "direct_light_";
    namestream << fSubRunNumbers[subRunIndex];
    WriteToRoot(&directLight,(namestream.str()).c_str());
    namestream.clear();
    namestream.str(std::string());

    namestream << "reflected_light_";
    namestream << fSubRunNumbers[subRunIndex];
    WriteToRoot(&reflectedLight,(namestream.str()).c_str());
    namestream.clear();
    namestream.str(std::string());

    namestream << "light_origin_";
    namestream << fSubRunNumbers[subRunIndex];
    WriteToRoot(&lightOrigin,(namestream.str()).c_str());
    namestream.clear();
    namestream.str(std::string());

    // Output
    detail << "DQTellieProc::FibreCheck:" << newline;
    detail << "- Direct light centre: x: " << directLight.X();
    detail << " y: " << directLight.Y();
    detail << " z: " << directLight.Z() << newline;
    detail << "- Reflected light centre: x: " << reflectedLight.X();
    detail << " y: " << reflectedLight.Y();
    detail << " z: " << reflectedLight.Z() << newline;
    detail  << "- Light Origin: x: " << lightOrigin.X();
    detail << " y: " << lightOrigin.Y();
    detail << " z: " << lightOrigin.Z() << newline;

    // Find closest fibre
    double minDist = 1e6;
    string likelyFibre = "";
    TVector3 likelyFibrePos(0,0,0);
    // There are 92 distinct fibre positions on PSUP. Fibres > 92 are spares and are in the same position as other fibres.
    // FIXME - loop over all patched fibres in RATDB instead (i.e. those connected to a TELLIE channel at time of run)
    int maxNumFibre = 95;
    for ( int i = 0; i <= maxNumFibre; i++)
    {
      stringstream fibNum;
      fibNum.str("");
      if    (i==46) fibNum << "FT047B";
      else if (i==95) fibNum << "FT101A";
      else      fibNum << "FT" << setw(3) << setfill('0') << i+1 << "A";
      DBLinkPtr fibre;
      TVector3 thisFibrePos(0,0,0);
      try { // Get fibre from DB
        fibre = DB::Get()->GetLink("FIBRE", fibNum.str());
        thisFibrePos.SetX(fibre->GetD("x"));
        thisFibrePos.SetY(fibre->GetD("y"));
        thisFibrePos.SetZ(fibre->GetD("z"));
      } catch ( DBNotFoundError &e ) {
        warn << "WARNING: DQTellieProc::FibreCheck " << fibNum.str();
        warn << " does not exist. Skipping." << newline;
        continue;
      }
      debug << "Fibre Pos: x: " << thisFibrePos.X();
      debug << " y: " << thisFibrePos.Y();
      debug << " z: " << thisFibrePos.Z() << newline;
      double dist = (lightOrigin-thisFibrePos).Mag();
      debug << "Dist from  fibre: "<<dist<<newline;
      if (dist < minDist)
      {
        minDist = dist;
        likelyFibrePos = thisFibrePos;
        likelyFibre = fibNum.str();
        debug << "Likely fibre: "<<likelyFibre<<newline;
      }
    }
    int dummy;
    TVector2 directProject = IcosProject(directLight, dummy);
    TVector2 reflectedProject = IcosProject(reflectedLight, dummy);
    TVector2 assumedProject = IcosProject(likelyFibrePos, dummy);
    TVector2 actualProject = IcosProject(fibrePos, dummy);

    // Draw 1-sigma contours around fitted light positions
    std::vector<TVector2> directContour;
    std::vector<TVector2> reflectedContour;
    DrawCircle(directLight, sigmaAngDir, directContour, 360);
    DrawCircle(reflectedLight, sigmaAngRef, reflectedContour, 360);
    PlotFibreCheck(directProject, reflectedProject, directContour, reflectedContour, assumedProject, actualProject, likelyFibre, subRunIndex);
    debug << "DQTellieProc::FibreCheck: Likely fibre: " << likelyFibre;
    debug << " distance from calculated origin " << minDist << newline;
    string firedFibre = fFibres[subRunIndex];

    // Check fibre fired is correct, account for 2 fibres at same node
    // FIXME - fibre names should not be hardcoded, instead should loop over nodes
    if (likelyFibre == firedFibre) {
      fCheckResult = true;
    } else if (likelyFibre == "FT022A" && firedFibre == "FT047B") {
      fCheckResult = true;
    } else if (likelyFibre == "FT047B" && firedFibre == "FT022A") {
      fCheckResult = true;
    } else if (likelyFibre == "FT026A" && firedFibre == "FT067A") {
      fCheckResult = true;
    } else if (likelyFibre == "FT067A" && firedFibre == "FT026A") {
      fCheckResult = true;
    } else if (likelyFibre == "FT091A" && firedFibre == "FT093A") {
      fCheckResult = true;
    } else if (likelyFibre == "FT093A" && firedFibre == "FT091A") {
      fCheckResult = true;
    } else {
      fCheckResult = false;
    }
    fFibreFiringGuess[subRunIndex] = likelyFibre;
    WriteToRoot(fTH2DPsupProjection[subRunIndex]);
    WriteToRoot(fTH2DPsupOffline[subRunIndex]);
    if(fCheckResult) {
      detail << "DQTellieProc::FibreCheck: run " << GetRunID();
      detail << " subrun " << fSubRunNumbers[subRunIndex];
      detail << " passed the fibre check (the correct fibre is firing)." << newline;
      fCorrectFibreCheck[subRunIndex] = 1;
    } else {
      info << "DQTellieProc::FibreCheck: run " << GetRunID();
      info << " subrun "<<fSubRunNumbers[subRunIndex];
      info << " failed the fibre check." << newline;
      info << "Fibre that should be firing: " << firedFibre;
      info << " - Fibre that appears to be firing: " << likelyFibre << newline;
      fCorrectFibreCheck[subRunIndex] = 0;
    }

  }

  /// Create detector view from above central point
  void DQTellieProc::FillHemisphere(const TVector3& center, std::vector<float> &occupancy, TGraph2D* graph) {

    const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();

    // Find rotation angles for this frame
    double rot_X, rot_Z;
    GetRotationAngles(center,rot_Z,rot_X);

    // Loop over all PMTs
    int counter=0;
    TVector3 pmtPos, newPos;
    for(unsigned int id=0; id<occupancy.size(); id++) {

      // Check PMTs are ok (FIXME - cuts should not be hardcoded)
      if(pmtInfo.GetType(id) != DU::PMTInfo::NORMAL) continue;
      pmtPos = pmtInfo.GetPosition(id);
      if (pmtPos.Mag()==0) continue;
      if (occupancy[id] < 0.001) continue; // less than 5 events (subrun)
      if (occupancy[id] > 0.1) continue; // more than 10% occupancy

      if (center.Angle(pmtPos) > TMath::Pi()/2.) continue;

      // Rotate particles into correct frame
      newPos = pmtPos;
      newPos.RotateZ(-rot_X);
      newPos.RotateY(-rot_Z);

      // Fill 2D graph
      if (newPos.Z() <= 0) continue; // not in hemisphere (safety check)
      double xpt = newPos.X()/1e3; // *cos(pi/4)/cos(newpos.Theta());
      double ypt = newPos.Y()/1e3; // *cos(pi/4)/cos(newpos.Theta());
      graph->SetPoint(counter, xpt, ypt, occupancy[id]);
      counter++;
    }

  }

  /// Fit 2D Gaussian surface over intensities for PMTs in a plane
  void DQTellieProc::FitLightSpot(TGraph2D* graph, double radius, double cone, double* params) {
    // Get input graph entries
    int npts = graph->GetN();
    double *xpts = graph->GetX();   // X projection
    double *ypts = graph->GetY();   // Y projection
    double *zpts = graph->GetZ();   // Intensity

    // First loop: Find maximum value in desired range
    double maxval=-1;
    for (int n=0; n<npts; n++) {
      double rpt = sqrt(xpts[n]*xpts[n] + ypts[n]*ypts[n]);
      if (rpt>radius) continue;
      double ang = asin(rpt/radius);
      if (ang/TMath::Pi()*180. > cone) continue;
      if (maxval<zpts[n]) maxval=zpts[n];
    }

    // Second loop: Fill values above 20% of maximum values into output graph
    TGraph2D *graf = new TGraph2D();
    int pts=0;
    for (int n=0; n<npts; n++) {
      if(zpts[n] < maxval*0.2) continue;
      double rpt = sqrt(xpts[n]*xpts[n] + ypts[n]*ypts[n]);
      if (rpt>radius) continue;
      double ang = asin(rpt/radius);
      if (ang/TMath::Pi()*180. > cone) continue;
      double scale = 1./cos(ang);
      graf->SetPoint(pts,scale*xpts[n],scale*ypts[n],zpts[n]);
      pts++;
    }

    // Fit 2D Gaussian surface to selected points
    TF2 *fit = new TF2("gaus2d","[0]*TMath::Gaus(x,[1],[3])*TMath::Gaus(y,[2],[3])",-10,10,-10,10);
    double aper = radius*tan(cone/2./180.*TMath::Pi()); // half input angle <=> 12 deg aperture
    fit->SetParameters(0.8*maxval,0.,0.,aper); // a priori: 80% max. intensity, centered, nominal aper
    fit->SetParLimits(0,0.2*maxval,maxval);  // amplitude range [20%, 100%] of max. intensity
    fit->SetParLimits(1,-aper,aper);       // x-pos range [-12, +12] degrees from weighted centre
    fit->SetParLimits(2,-aper,aper);       // y-pos range [-12, +12] degrees from weighted centre
    fit->SetParLimits(3,0.5*aper,1.5*aper);  // sigma range [-50%, +50%] of nominal aperture
    graf->Fit("gaus2d","R,q");         // fit gaussian (force range, quiet mode)

    // Raw fit parameters
    double amp = fit->GetParameter(0);
    double mux = fit->GetParameter(1);
    double muy = fit->GetParameter(2);
    double sig = fit->GetParameter(3);

    // Scale mean value back to sphere
    double rad = sqrt(mux*mux+muy*muy);
    double sx = mux*cos(atan(rad/radius));
    double sy = muy*cos(atan(rad/radius));

    // Scale (mean +- sigma) values back to sphere
    double xp = (mux+sig);
    double rxp = sqrt(xp*xp+muy*muy);
    double sxp = xp*cos(atan(rxp/radius));
    double xm = (mux-sig);
    double rxm = sqrt(xm*xm+muy*muy);
    double sxm = xm*cos(atan(rxm/radius));
    double yp = (muy+sig);
    double ryp = sqrt(yp*yp+mux*mux);
    double syp = yp*cos(atan(ryp/radius));
    double ym = (muy-sig);
    double rym = sqrt(ym*ym+mux*mux);
    double sym = ym*cos(atan(rym/radius));

    // Average over all scaled sigmas
    double sigma = (fabs(sxp-sx)+fabs(sxm-sx)+fabs(syp-sy)+fabs(sym-sy))/4.;
    //printf("Fit result scaled to sphere:\tp1 = %.3lf\tp2 = %.3lf\tp3 = %.3lf\n",sx,sy,sigma);

    // Output fit results
    params[0] = amp;  // amplitude
    params[1] = sx;   // mu_x
    params[2] = sy;   // mu_y
    params[3] = sigma;  // sigma

    if (graf) delete graf;
    if (fit) delete fit;
  }

  // Find out if direct light from fibre is affected by belly plates
  bool DQTellieProc::CheckBellyPlate(int subRunIndex) {

    // The fibres affected by belly plates can be hardcoded for now
    std::vector<string> belly_fibres;
    belly_fibres.push_back("FT059A");
    belly_fibres.push_back("FT066A");
    belly_fibres.push_back("FT065A");
    belly_fibres.push_back("FT069A");
    belly_fibres.push_back("FT068A");
    belly_fibres.push_back("FT074A");
    belly_fibres.push_back("FT073A");
    belly_fibres.push_back("FT038A");
    belly_fibres.push_back("FT039A");
    belly_fibres.push_back("FT032A");
    belly_fibres.push_back("FT033A");
    belly_fibres.push_back("FT042A");
    belly_fibres.push_back("FT041A");
    belly_fibres.push_back("FT048A");
    belly_fibres.push_back("FT101A");
    belly_fibres.push_back("FT051A");
    belly_fibres.push_back("FT050A");
    belly_fibres.push_back("FT057A");
    belly_fibres.push_back("FT056A");
    belly_fibres.push_back("FT060A");

    // Check if fibre firing for this subrun is affected by belly plates
    string fibre = fFibres[subRunIndex];
    if (find(belly_fibres.begin(), belly_fibres.end(), fibre) != belly_fibres.end()) {
      return true;
    }
    return false;
  }

  /// Draw a circle around a point in a plane orthogonal to an angle with N dots
  void DQTellieProc::DrawCircle(const TVector3& center, double angle,
                  std::vector<TVector2>& dots, int NDOTS) {

    int dummy;
    TVector3 dot(0,0,0);
    for (int i=0; i<NDOTS; i++) {
      // Draw circle around center point
      dot.SetMagThetaPhi(center.Mag(), TMath::Pi()*angle/180., i*2.*TMath::Pi()/NDOTS);
      dot.RotateUz(center.Unit());   // rotate into frame of input vector

      // Store projected circle in vector
      dots.push_back(IcosProject(dot,dummy));
    }
  }

  /// Get rotation angles to rotate from given direction to xyz-frame
  void DQTellieProc::GetRotationAngles(const TVector3& vec, double &rot_Z, double &rot_X) {

    // Unit vectors
    const TVector3 e1(1,0,0);
    const TVector3 e2(0,1,0);
    const TVector3 e3(0,0,1);

    // Get rotation angle from z-axis to central direction
    rot_Z = acos(e3*(vec.Unit()));     // rotate counter-clockwise, towards x-axis by [0,pi]

    // Get rotation angle from x-axis to central direction (projected to transverse plane)
    int sign = 0;
    if (vec[1] > 0) sign = +1;
    else      sign = -1;
    TVector3 vect(vec[0],vec[1],0);
    rot_X = sign*acos(e1*(vect.Unit()));   // rotate around z-axis by [-pi,pi]
  }

  /// Display vector as a string
  string DQTellieProc::printVector(const TVector3& v) {
    string out;
    if (v.Mag() < 10) out = Form("(%.3f, %.3f, %.3f)", v.X(),  v.Y(), v.Z());
    else        out = Form("(%.1f, %.1f, %.1f)", v.X(),  v.Y(), v.Z());
    return out.c_str();
  }

  void DQTellieProc::PlotTac(std::vector<float> xPeaks,int subRunIndex)
  {
    gStyle->SetOptStat(0);
    TCanvas *C1 = new TCanvas("C1");
    C1->cd()->SetLogy();

    // Draw TAC distributions
    float maxval = fTH1DTacCount[subRunIndex]->GetMaximum();
    fTH1DTacCount[subRunIndex]->Draw();
    fTH1DTacCount[subRunIndex]->SetXTitle("TAC");
    fTH1DTacCount[subRunIndex]->SetYTitle("Counts");
    fTH1DTacCount[subRunIndex]->GetXaxis()->SetTitleOffset(1.4);
    fTH1DTacCount[subRunIndex]->GetYaxis()->SetTitleOffset(1.4);
    fTH1DTacCount[subRunIndex]->GetYaxis()->SetRangeUser(maxval/2e3, maxval*2);
    fTH1DTacCount[subRunIndex]->SetLineColor(1);
    fTH1DTacCountDirect[subRunIndex]->Draw("same");
    fTH1DTacCountDirect[subRunIndex]->SetLineColor(2);
    fTH1DTacCountReflected[subRunIndex]->Draw("same");
    fTH1DTacCountReflected[subRunIndex]->SetLineColor(4);

    // Draw peaks found
    TMarker marker;
    marker.SetMarkerStyle(22);
    marker.SetMarkerColor(2);
    marker.SetMarkerSize(1.5);
    for (vector<float>::iterator it = xPeaks.begin(); it != xPeaks.end(); it++) {
      float xp = *it;
      int bin = fTH1DTacCount[subRunIndex]->GetXaxis()->FindBin(xp);
      float yp = fTH1DTacCount[subRunIndex]->GetBinContent(bin);
      marker.DrawMarker(xp,yp);
    }

    // Draw legend
    float x0=0.13, y0=0.70, x1=0.43, y1=0.88;
    int maxbin = fTH1DTacCount[subRunIndex]->GetMaximumBin();
    float xmax = fTH1DTacCount[subRunIndex]->GetXaxis()->GetBinCenter(maxbin);
    if (xmax <= 2048) { // move legend to top right
      x0=0.57; x1=0.87;
    }
    TLegend leg(x0,y0,x1,y1);
    string title = Form("Total (%s) ", fFibres[subRunIndex].c_str());
    leg.AddEntry(fTH1DTacCount[subRunIndex],title.c_str(),"l");
    leg.AddEntry(fTH1DTacCountDirect[subRunIndex],"Direct light (>165#circ) ","l");
    leg.AddEntry(fTH1DTacCountReflected[subRunIndex],"Reflected light (<15#circ) ","l");
    leg.SetBorderSize(0);
    leg.Draw();

    // Save graph
    gPad->SetTicks(1,1);
    C1->Update();
    string filename = (fTH1DTacCount[subRunIndex])->GetName();
    filename = filename+".png";
    stringstream namestream;
    namestream << "TAC Peaks subrun " << subRunIndex;
    C1->SetName(namestream.str().c_str());
    WriteToRoot(C1);
    C1->Print(filename.c_str());
    delete C1;
    C1 = NULL;
  }

  void DQTellieProc::GetTopIndices(const std::vector<float> x, const std::vector<float> y,
                   size_t& early, size_t& peak, size_t& late) {
    if(x.size()==0 || y.size()==0) return;

    // Sort vector and find highest three peaks
    std::vector<float> w = y;
    std::sort(w.begin(), w.end());
    float first  = *(w.end()-1); // highest peak
    float second = *(w.end()-2); // second highest peak
    float third  = *(w.end()-3); // third highest peak

    // Find indices for highest three peaks
    size_t idx = 0;
    for(vector<float>::const_iterator it = y.begin(); it != y.end(); it++) {
      if (*it == first) peak = idx;
      if (*it == second) early = idx;
      if (*it == third) late = idx;
      idx++;
    }
  }

  void DQTellieProc::TimingChecks(int subRunIndex)
  {
    // Estimate linear background
    TF1 *line = new TF1("line","pol1",0,4025);
    fTH1DTacCountReflected[subRunIndex]->Fit("line","qn"); //qn = quiet & no draw

    // Find peaks in reflected/direct light spectra
    TSpectrum *peakFinderRef = new TSpectrum(2*fNPeaks[subRunIndex]);
    int nFoundRef = peakFinderRef->Search(fTH1DTacCountReflected[subRunIndex], fTACReflectedSigmaThresh, "new", fTACReflectedHeightThresh);
    TSpectrum *peakFinderDir = new TSpectrum(2*fNPeaks[subRunIndex]);
    int nFoundDir = peakFinderDir->Search(fTH1DTacCountDirect[subRunIndex], fTACPromptSigmaThresh, "new", fTACPromptHeightThresh);
    int nFound = nFoundDir + nFoundRef;
    debug << "Number of peaks found: " << nFound << newline;
    vector<double> par(nFound*3+2);
    par[0] = line->GetParameter(0);
    par[1] = line->GetParameter(1);
#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
    double* xPeaks_checkDir = peakFinderDir->GetPositionX();
    double* xPeaks_checkRef = peakFinderRef->GetPositionX();
#else
    float* xPeaks_checkDir = peakFinderDir->GetPositionX();
    float* xPeaks_checkRef = peakFinderRef->GetPositionX();
#endif
    vector<float> xPeaks;
    vector<float> yPeaks;
    fNPeaks[subRunIndex] = 0;
    debug  << "Checking peaks" << newline;
    for(int i=0; i < nFoundDir; i++) {
      float xp = xPeaks_checkDir[i];
      int bin = fTH1DTacCountDirect[subRunIndex]->GetXaxis()->FindBin(xp);
      float yp = fTH1DTacCountDirect[subRunIndex]->GetBinContent(bin);
      xPeaks.push_back(xp);
      yPeaks.push_back(yp);
      par[3*fNPeaks[subRunIndex]+2] = yp;
      par[3*fNPeaks[subRunIndex]+3] = xp;
      par[3*fNPeaks[subRunIndex]+4] = 200.;
      fNPeaks[subRunIndex]++;
      debug << "Direct y:" << yp << "x:" << xp << newline;
    }
    for(int i=0; i < nFoundRef; i++) {
      float xp = xPeaks_checkRef[i];
      int bin = fTH1DTacCountReflected[subRunIndex]->GetXaxis()->FindBin(xp);
      float yp = fTH1DTacCountReflected[subRunIndex]->GetBinContent(bin);
      xPeaks.push_back(xp);
      yPeaks.push_back(yp);
      par[3*fNPeaks[subRunIndex]+2] = yp;
      par[3*fNPeaks[subRunIndex]+3] = xp;
      par[3*fNPeaks[subRunIndex]+4] = 200.;
      fNPeaks[subRunIndex]++;
      debug << "Reflected y:" << yp << "x:" << xp << newline;
    }

    // Find peak order
    size_t pre_idx = -1;
    size_t prompt_idx = -1;
    size_t late_idx = -1;
    GetTopIndices(xPeaks,yPeaks,pre_idx,prompt_idx,late_idx);
    /*
    MinMaxIndex(xPeaks,late_idx,pre_idx);
    size_t prompt_idx = 3 - pre_idx - late_idx;
    */
    PlotTac(xPeaks,subRunIndex);

    //Check 7 - Number of timing peaks found
    if (fNPeaks[subRunIndex] == 3) {
      detail << "DQTellieProc::TimingChecks: run " << GetRunID();
      detail << " subrun " << fSubRunNumbers[subRunIndex];
      detail << " passed the number of timing peaks check" << newline;
      fCheckResult = true;
    } else {
      info << "DQTellieProc::TimingChecks: run " << GetRunID();
      info << " subrun " << fSubRunNumbers[subRunIndex];
      info << " failed the number of timing peaks check";
      info << " (3 peaks expected, but " << fNPeaks[subRunIndex];
      info << " found.)" << newline;
      fCheckResult = false;
    }
    if(fCheckResult) fPeakNumberCheck[subRunIndex] = 1;
    else       fPeakNumberCheck[subRunIndex] = 0;

    //Check 8 - Peak amplitude check
    std::stringstream namestream;
    namestream << "amp_pre_subrun_";
    namestream << fSubRunNumbers[subRunIndex];
    TVectorD amp_pre(1);
    amp_pre[0] = par[3*pre_idx+2];
    WriteToRoot(&amp_pre, (namestream.str()).c_str() );
    namestream.clear();
    namestream.str(std::string());

    namestream << "amp_prompt_subrun_";
    namestream << fSubRunNumbers[subRunIndex];
    TVectorD amp_prompt(1);
    amp_prompt[0] = par[3*prompt_idx+2];
    WriteToRoot(&amp_prompt, (namestream.str()).c_str() );
    namestream.clear();
    namestream.str(std::string());

    namestream << "amp_late_subrun_";
    namestream << fSubRunNumbers[subRunIndex];
    TVectorD amp_late(1);
    amp_late[0] = par[3*late_idx+2];
    WriteToRoot(&amp_late, (namestream.str()).c_str() );
    namestream.clear();
    namestream.str(std::string());

    if (fTACPromptHeightRatio*amp_pre[0] < amp_prompt[0] &&
      fTACPromptHeightRatio*amp_late[0] < amp_prompt[0]) {
      detail << "DQTellieProc::TimingChecks: run " << GetRunID();
      detail << " subrun " << fSubRunNumbers[subRunIndex];
      detail << " passed the timing peak amplitude check." << newline;
      fCheckResult = true;
    } else {
      info << "DQTellieProc::TimingChecks: run " << GetRunID();
      info << " subrun " << fSubRunNumbers[subRunIndex];
      info << " failed the timing peak amplitude check";
      info << " (Prompt peak must be at least " << fTACPromptHeightRatio;
      info << " times bigger)." << newline;
      info << " Prompt amplitude: " << amp_prompt[0];
      info << " pre pulse amplitude " << amp_pre[0];
      info << " late pulse amplitude " << amp_late[0] << newline;
      fCheckResult = false;
    }
    if(fCheckResult) fPromptPeakAmplitudeCheck[subRunIndex] = 1;
    else             fPromptPeakAmplitudeCheck[subRunIndex] = 0;

    // Check 9 - Relative peak time check
    namestream << "tac_prompt_subrun_";
    namestream << fSubRunNumbers[subRunIndex];
    TVectorD tac_prompt(1);
    tac_prompt[0]  = par[3*prompt_idx+3];
    WriteToRoot(&tac_prompt, (namestream.str()).c_str() );
    namestream.clear();
    namestream.str(std::string());

    if (tac_prompt[0] >= fTACPromptPeakLowerLimit &&
      tac_prompt[0] <= fTACPromptPeakUpperLimit) {
      detail << "DQTellieProc::TimingChecks: run " << GetRunID();
      detail << " subrun "<< fSubRunNumbers[subRunIndex];
      detail << " passed the prompt peak time check." << newline;
      fCheckResult = true;
    } else {
      info << "DQTellieProc::TimingChecks: run " << GetRunID();
      info << " subrun " << fSubRunNumbers[subRunIndex];
      info << " failed the prompt peak time check";
      info << " (value not between " << fTACPromptPeakLowerLimit;
      info << " and " << fTACPromptPeakUpperLimit << " ADC)." << newline;
      info << "Prompt peak time " << tac_prompt[0] << " ADC." << newline;
      fCheckResult = false;
    }
    if(fCheckResult) fPromptPeakADCCountCheck[subRunIndex] = 1;
    else             fPromptPeakADCCountCheck[subRunIndex] = 0;

    // Check 10 - Relative peak time check
    namestream << "tac_late_subrun_";
    namestream << fSubRunNumbers[subRunIndex];
    TVectorD tac_late(1);
    tac_late[0] = par[3*late_idx+3];
    WriteToRoot(&tac_late, (namestream.str()).c_str() );
    namestream.clear();
    namestream.str(std::string());

    namestream << "tac_pre_subrun_";
    namestream << fSubRunNumbers[subRunIndex];
    TVectorD tac_pre(1);
    tac_pre[0]  = par[3*pre_idx+3];
    WriteToRoot(&tac_pre, (namestream.str()).c_str() );
    namestream.clear();
    namestream.str(std::string());

    if (tac_late[0] >= (tac_prompt[0]-fTACLateOffsetUpper) &&
      tac_late[0] <= (tac_prompt[0]-fTACLateOffsetLower) &&
      tac_pre[0] >= (tac_prompt[0]+fTACEarlyOffsetLower) &&
      tac_pre[0] <= (tac_prompt[0]+fTACEarlyOffsetUpper)) {
      detail << "DQTellieProc::TimingChecks: run " << GetRunID();
      detail << " subrun "<< fSubRunNumbers[subRunIndex];
      detail << " passed the peak time check." << newline;
      fCheckResult = true;
    } else {
      info << "DQTellieProc::TimingChecks: run " << GetRunID();
      info << " subrun "<< fSubRunNumbers[subRunIndex];
      info << " failed the peak time check";
      info << " (pre pulse not between " << fTACEarlyOffsetLower;
      info << " and " << fTACEarlyOffsetUpper;
      info << " ADC more than prompt; or";
      info << " late pulse not between " << fTACLateOffsetLower;
      info << " and " << fTACLateOffsetUpper;
      info << " ADC less than prompt.)" << newline;
      info << "Prompt time: " << tac_prompt[0] << " ADC.";
      info << " Pre pulse time " << tac_pre[0] << " ADC.";
      info << " Late pulse time " << tac_late[0] << " ADC." << newline;
      fCheckResult = false;
    }
    if(fCheckResult) fPeakTimeSpacingCheck[subRunIndex] = 1;
    else             fPeakTimeSpacingCheck[subRunIndex] = 0;

    /*
    // Output for tuning cut parameters only
    info << "PEAKSTATS:";
    info << " " << tac_prompt[0] << " " << amp_prompt[0];
    info << " " << tac_pre[0] << " " << amp_pre[0];
    info << " " << tac_late[0] << " " << amp_late[0];
    info << newline;
    */

    fPromptPeakADCCount[subRunIndex] = tac_prompt[0];
    fPrePeakADCCount[subRunIndex] = tac_pre[0];
    fLatePeakADCCount[subRunIndex] = tac_late[0];
    WriteToRoot(fTH1DTacCount[subRunIndex]);
    WriteToRoot(fTH1DTacCountDirect[subRunIndex]);
    WriteToRoot(fTH1DTacCountReflected[subRunIndex]);
    WriteToRoot(fTH1DCalibTime[subRunIndex]);
    WriteToRoot(fTH1DPulseSeparation[subRunIndex]);
    WriteToRoot(fTH1DPulseSeparationStolen[subRunIndex]);

    TCanvas *C1 = NULL;
    string filename;
    double minval, maxval;

    // Plot calibrated time
    C1 = new TCanvas("C1");
    fTH1DCalibTime[subRunIndex]->SetXTitle("Time (ns)");
    fTH1DCalibTime[subRunIndex]->SetYTitle("Counts");
    fTH1DCalibTime[subRunIndex]->GetXaxis()->SetTitleOffset(1.4);
    fTH1DCalibTime[subRunIndex]->GetYaxis()->SetTitleOffset(1.4);
    filename = fTH1DCalibTime[subRunIndex]->GetName();
    filename = filename+".png";
    C1->cd()->SetLogy();
    fTH1DCalibTime[subRunIndex]->Draw();
    fTH1DCalibTime[subRunIndex]->Write();
    minval = fTH1DCalibTime[subRunIndex]->GetBinCenter(
             fTH1DCalibTime[subRunIndex]->FindFirstBinAbove(0) - 10 );
    maxval = fTH1DCalibTime[subRunIndex]->GetBinCenter(
             fTH1DCalibTime[subRunIndex]->FindLastBinAbove(0) + 10 );
    fTH1DCalibTime[subRunIndex]->GetXaxis()->SetRangeUser(minval,maxval);
    gPad->SetTicks(1,1);
    C1->Update();
    C1->Print( filename.c_str() );
    delete C1;
    C1 = NULL;

    // Plot pulse separations
    C1 = new TCanvas("C1");
    fTH1DPulseSeparation[subRunIndex]->SetXTitle("TELLIE event separation (ms)");
    fTH1DPulseSeparation[subRunIndex]->SetYTitle("Counts");
    fTH1DPulseSeparation[subRunIndex]->GetXaxis()->SetTitleOffset(1.4);
    fTH1DPulseSeparation[subRunIndex]->GetYaxis()->SetTitleOffset(1.4);
    filename = fTH1DPulseSeparation[subRunIndex]->GetName();
    filename = filename+".png";
    C1->cd()->SetLogy();
    fTH1DPulseSeparation[subRunIndex]->SetLineColor(4);
    fTH1DPulseSeparation[subRunIndex]->SetFillStyle(3145);
    fTH1DPulseSeparation[subRunIndex]->SetFillColor(4);
    fTH1DPulseSeparation[subRunIndex]->Draw();
    fTH1DPulseSeparation[subRunIndex]->Write();
    fTH1DPulseSeparationStolen[subRunIndex]->SetLineColor(2);
    fTH1DPulseSeparationStolen[subRunIndex]->SetFillStyle(3154);
    fTH1DPulseSeparationStolen[subRunIndex]->SetFillColor(2);
    fTH1DPulseSeparationStolen[subRunIndex]->Draw("same");
    fTH1DPulseSeparationStolen[subRunIndex]->Write();
    minval = fTH1DPulseSeparation[subRunIndex]->GetBinCenter(
             fTH1DPulseSeparation[subRunIndex]->FindFirstBinAbove(0) - 10 );
    maxval = fTH1DPulseSeparation[subRunIndex]->GetBinCenter(
             fTH1DPulseSeparation[subRunIndex]->FindLastBinAbove(0) + 10 );
    fTH1DPulseSeparation[subRunIndex]->GetXaxis()->SetRangeUser(minval,maxval);
    fTH1DPulseSeparation[subRunIndex]->GetYaxis()->SetRangeUser(0.5,1e4);
    TLegend *leg = new TLegend(0.13, 0.78, 0.43, 0.86);
    leg->AddEntry(fTH1DPulseSeparation[subRunIndex],"EXTA triggers","f");
    leg->AddEntry(fTH1DPulseSeparationStolen[subRunIndex],"Stolen triggers","f");
    leg->SetBorderSize(0);
    leg->Draw();
    gStyle->SetOptStat(101111); // include overflow bins
    gPad->SetTicks(1,1);
    C1->Update();
    C1->Print( filename.c_str() );
    delete C1;
    C1 = NULL;

    info << "Finished Timing Checks" << newline;
  }

  void DQTellieProc::EndOfRun( DS::Run& run )
  {
    //Doing checks for each subrun
    for(unsigned int s=0; s<fSubRunNumbers.size(); s++) {

      // Check 1 - trigger check
      double trigCheckResult = 100.0*(static_cast<double>(fTellieCount[s])/
                      static_cast<double>(fTellieEvents[s]));
      detail << "tellie count for subrun " << fSubRunNumbers[s];
      detail << ":  " << fTellieCount[s];
      detail << " tellie events: " << fTellieEvents[s] << newline;
      if( trigCheckResult >= fTrigCheckThresh && trigCheckResult <= 100.0) {
        detail << "DQTellieProc::EndOfRun: run " << GetRunID();
        detail << " subrun " << fSubRunNumbers[s];
        detail << " passed trigger check" << newline;
        fCheckResult = true;
      } else if( trigCheckResult > 100.0 ) {
        info << "DQTellieProc::EndOfRun: run " << GetRunID();
        info << " subrun " << fSubRunNumbers[s];
        info << " failed trigger check" << newline;
        info << " --> " << trigCheckResult;
        info << "% of events had the correct active triggers.";
        info << " Number of tellie events in check is wrong!" << newline;
        fCheckResult = false;
      } else if ( fTellieEvents[s] - fTellieCount[s] == fStolenTriggerCount[s] ) {
        detail << "DQTellieProc::EndOfRun: run " << GetRunID();
        detail << " subrun " << fSubRunNumbers[s];
        detail << " passed trigger check (";
        detail << fStolenTriggerCount[s] << " events stolen by other triggers were recovered)." << newline;
        fCheckResult = true;
      } else {
        info << "DQRunProc::EndOfRun: run " << GetRunID();
        info << " subrun " << fSubRunNumbers[s];
        info << " failed trigger check" << newline;
        info << " --> " << trigCheckResult;
        info << "% of events (less than " << fTrigCheckThresh;
        info << "% threshold) had the correct active triggers" << newline;
        info << (fTellieEvents[s] - fTellieCount[s]);
        info << " TELLIE events did not trigger on EXTA (";
        info << fStolenTriggerCount[s] << " events stolen by other triggers were recovered)." << newline;
        warn << "WARNING - " << fTellieEvents[s] - fTellieCount[s] - fStolenTriggerCount[s];
        warn << " TELLIE events went missing!" << newline;
        fCheckResult = false;
      }
      if (fCheckResult) fTriggerCheck[s] = 1;
      else              fTriggerCheck[s] = 0;

      // Check 2 - run time checks
      double runLength = (fLastEventTime[s] - fFirstEventTime[s]).GetSeconds() +
                 ((fLastEventTime[s] - fFirstEventTime[s]).GetNanoSeconds()*1e-9);
      debug << "DQTellieProc::EndOfRun: for run " << GetRunID();
      debug << " subrun " << fSubRunNumbers[s] << " time from 10 Mhz Clock";
      debug << " (last (= " << fLastEventTime[s].GetSeconds() << ") -";
      debug << " first (= " << fFirstEventTime[s].GetSeconds() << "))";
      debug << " = " << runLength << newline;
      fSubRunLengths[s] = runLength;
      double minRunLength = fTellieEvents[s]/fTelliePulseRate[s];

      if( runLength > minRunLength ) {
        detail << "DQTellieProc::EndOfRun: run " << GetRunID();
        detail << " subrun " << fSubRunNumbers[s];
        detail << " (run time = " << ( runLength );
        detail << ") passed max/min run time check" << newline;
        fCheckResult = true;
      } else {
        info << "DQTellieProc::EndOfRun: run " << GetRunID();
        info << " subrun " << fSubRunNumbers[s];
        info << " failed min run length check" << newline;
        info << " --> run length (" << ( runLength );
        info << ") was less than minimum run length (";
        info << minRunLength << ")" << newline;
        fCheckResult = false;
      }
      if (fCheckResult) fSubRunLengthsCheck[s] = 1;
      else              fSubRunLengthsCheck[s] = 0;

      //Check 3 - pulse separation checks
      double pulseSeparationCheck = 100.*(1.0 - static_cast<double>(fFailPulseSeparationCount[s]) /
                                           static_cast<double>(fTellieEvents[s]) );
      fPulseSeparationEfficiency[s] = pulseSeparationCheck;
      if (pulseSeparationCheck >= fPulseSeparationThresh) {
        detail << "DQTellieProc::EndOfRun: run " << GetRunID();
        detail <<" subrun " << fSubRunNumbers[s];
        detail << " (pulse separation efficiency = " << pulseSeparationCheck;
        detail << "%) passed  max/min pulse separation check" << newline;
        fCheckResult = true;
      } else {
        info << "DQTellieProc::EndOfRun: run " << GetRunID();
        info << " subrun " << fSubRunNumbers[s];
        info << " failed pulse separation check";
        info << " --> pulse separation efficiency (" << pulseSeparationCheck;
        info << "%) was less than minimum efficiency for min/max pulse separation (";
        info << fPulseSeparationThresh << "%)" << newline;
        fCheckResult = false;
      }
      if (fCheckResult) fPulseSeparationEfficiencyCheck[s] = 1;
      else              fPulseSeparationEfficiencyCheck[s] = 0;

      //Check 4 - NHit
      if ( fFailMaxNHitCount[s] == 0) {
        detail << "DQTellieProc::EndOfRun: run " << GetRunID();
        detail << " subrun "<<fSubRunNumbers[s];
        detail << " passed the maximum NHit check" << newline;
        fCheckResult = true;
      } else {
        info << "DQTellieProc::EndOfRun: run" << GetRunID();
        info << " subrun "<<fSubRunNumbers[s];
        info << " failed NHit check --> " << fFailMaxNHitCount[s];
        info << " events had more than the maximum NHit (";
        info << fMaxNHit << ")." << newline;
        fCheckResult = false;
      }
      if(fCheckResult) fMaxNHitCheck[s] = 1;
      else             fMaxNHitCheck[s] = 0;

      //Check 5 - Average NHit
      //Doing this to avoid nan being output into RATDB file
      //If there are no tellie events then the average NHit for those events is 0
      if(fTellieCount[s] == 0) fAvgNHit[s] = 0;
      else           fAvgNHit[s] /= static_cast<double>(fTellieCount[s]);
      if ( fAvgNHit[s] < fMaxAvgNHit) {
        detail << "DQTellieProc::EndOfRun: run " << GetRunID();
        detail << " subrun " << fSubRunNumbers[s];
        detail << " (average NHit = " << fAvgNHit[s];
        detail << ") passed the maximum average NHit check" << newline;
        fCheckResult = true;
      } else {
        info << "DQTellieProc::EndOfRun: run " << GetRunID();
        info << " subrun "<< fSubRunNumbers[s];
        info << " failed average NHit check";
        info << " --> Average NHit (" << fAvgNHit[s];
        info << ") was more than the maximum average NHit (";
        info << fMaxAvgNHit << ")." << newline;
        fCheckResult = false;
      }
      if(fCheckResult) fAvgNHitCheck[s] = 1;
      else             fAvgNHitCheck[s] = 0;

      // Check 6 - Fibre
      info << "Doing Fibre  check for subrun: "<< fSubRunNumbers[s] << newline;
      FibreCheck(s);

      // Checks 7-10 - Timing checks
      info << "Doing Timing check for subrun: "<< fSubRunNumbers[s] << newline;
      TimingChecks(s);

      // Save the NHit Histos
      WriteToRoot(fTH1INHits[s]);
      TCanvas *C1 = new TCanvas("C1");
      C1->cd();
      stringstream namestream;
      namestream << "NHit Distribution "<< fSubRunNumbers[s];
      C1->SetName(namestream.str().c_str());
      fTH1INHits[s]->Draw();
      int lastbin = fTH1INHits[s]->FindLastBinAbove(0);
      int maxval = fTH1INHits[s]->GetXaxis()->GetBinUpEdge(lastbin);
      fTH1INHits[s]->GetXaxis()->SetRangeUser(0,50*TMath::Ceil(maxval/50.));
      fTH1INHits[s]->GetXaxis()->SetTitle("NHits");
      fTH1INHits[s]->GetYaxis()->SetTitle("Counts");
      fTH1INHits[s]->GetXaxis()->SetTitleOffset(1.4);
      fTH1INHits[s]->GetYaxis()->SetTitleOffset(1.4);
      gPad->SetTicks(1,1);
      C1->Update();
      string filename = (fTH1INHits[s])->GetName();
      filename = filename+".png";
      C1->Print(filename.c_str());
      delete C1;
      C1 = NULL;

    }

    //Adding the vectors to RATDB
    AddToRATDB("subrun_numbers",fSubRunNumbers);
    AddToRATDB("subrun_run_times",fSubRunLengths);
    AddToRATDB("expected_tellie_events",fTellieEvents);
    AddToRATDB("actual_tellie_events",fTellieCount);
    AddToRATDB("pulse_separation_efficiency",fPulseSeparationEfficiencyCheck);
    AddToRATDB("peak_numbers",fNPeaks);
    AddToRATDB("prompt_peak_adc_count",fPromptPeakADCCount);
    AddToRATDB("pre_peak_adc_count",fPrePeakADCCount);
    AddToRATDB("late_peak_adc_count",fLatePeakADCCount);
    AddToRATDB("average_nhit",fAvgNHit);
    AddToRATDB("more_max_nhit_events",fFailMaxNHitCount);
    AddToRATDB("run_mode",fRunMode);
    AddToRATDB("fibre_firing",fFibres);
    AddToRATDB("fibre_firing_guess",fFibreFiringGuess);

    //Add Check vectors to RATDB
    AddToRATDB("peak_number_check",fPeakNumberCheck);
    AddToRATDB("prompt_peak_adc_count_check",fPromptPeakADCCountCheck);
    AddToRATDB("prompt_peak_amplitude_check",fPromptPeakAmplitudeCheck);
    AddToRATDB("adc_peak_time_spacing_check",fPeakTimeSpacingCheck);
    AddToRATDB("trigger_check",fTriggerCheck);
    AddToRATDB("subrun_run_length_check",fSubRunLengthsCheck);
    AddToRATDB("pulse_separation_efficiency_check",fPulseSeparationEfficiencyCheck);
    AddToRATDB("max_nhit_check",fMaxNHitCheck);
    AddToRATDB("avg_nhit_check",fAvgNHitCheck);
    AddToRATDB("correct_fibre_check",fCorrectFibreCheck);

    //Check if all subruns pass if so pass run, otherwize fail
    bool peakNumberRunCheck = true;
    bool promptPeakADCRunCheck = true;
    bool promptPeakAmplitudeRunCheck = true;
    bool adcPeakTimeSpacingRunCheck = true;
    bool triggerRunCheck = true;
    bool subrunLengthsRunCheck = true;
    bool pulseSeparationEfficiencyRunCheck = true;
    bool maxNHitRunCheck = true;
    bool avgNHitRunCheck = true;
    bool correctFibreRunCheck = true;

    for(unsigned int s=0; s<fSubRunNumbers.size(); s++){
      if(peakNumberRunCheck){
        if(fPeakNumberCheck[s] == 0){
          peakNumberRunCheck = false;
        }
      }

      if(promptPeakADCRunCheck){
        if(fPromptPeakADCCountCheck[s] == 0){
          promptPeakADCRunCheck  = false;
        }
      }

      if(promptPeakAmplitudeRunCheck){
        if(fPromptPeakAmplitudeCheck[s] == 0){
          promptPeakAmplitudeRunCheck  = false;
        }
      }

      if(adcPeakTimeSpacingRunCheck){
        if(fPeakTimeSpacingCheck[s] == 0){
          adcPeakTimeSpacingRunCheck  = false;
        }
      }

      if(triggerRunCheck){
        if(fTriggerCheck[s] == 0){
          triggerRunCheck = false;
        }
      }

      if(subrunLengthsRunCheck){
        if(fSubRunLengthsCheck[s] == 0){
          subrunLengthsRunCheck = false;
        }
      }

      if(pulseSeparationEfficiencyRunCheck){
        if(fPulseSeparationEfficiencyCheck[s] == 0){
          pulseSeparationEfficiencyRunCheck = false;
        }
      }

      if(maxNHitRunCheck){
        if(fMaxNHitCheck[s] == 0){
          maxNHitRunCheck = false;
        }
      }

      if(avgNHitRunCheck){
        if(fAvgNHitCheck[s] == 0){
          avgNHitRunCheck = false;
        }
      }

      if(correctFibreRunCheck){
        if(fCorrectFibreCheck[s] == 0){
          correctFibreRunCheck = false;
        }
      }
    }

    Update("peak_number", peakNumberRunCheck);
    Update("peak_amplitude", promptPeakAmplitudeRunCheck);
    Update("prompt_time", promptPeakAmplitudeRunCheck);
    Update("peak_time", adcPeakTimeSpacingRunCheck);
    Update( "trigger", triggerRunCheck );
    Update( "run_length", subrunLengthsRunCheck );
    Update( "pulse_separation", pulseSeparationEfficiencyRunCheck);
    Update( "max_nhit", maxNHitRunCheck);
    Update( "avg_nhit", avgNHitRunCheck);
    Update("fibre",correctFibreRunCheck);
    for(unsigned int s=0; s<fSubRunNumbers.size(); s++){
      delete fTH1DLcnHits[s];
      fTH1DLcnHits[s] = NULL;
      delete fTH1DTacCount[s];
      fTH1DTacCount[s] = NULL;
      delete fTH1DCalibTime[s];
      fTH1DCalibTime[s] = NULL;
      delete fTH1DPulseSeparation[s];
      fTH1DPulseSeparation[s] = NULL;
      delete fTH1DPulseSeparationStolen[s];
      fTH1DPulseSeparationStolen[s] = NULL;
      delete fTH2DPsupProjection[s];
      fTH2DPsupProjection[s] = NULL;
      delete fTH2DPsupOffline[s];
      fTH2DPsupOffline[s] = NULL;
      delete fTH1INHits[s];
      fTH1INHits[s] = NULL;
    }
    DataQualityProc::EndOfRun( run );
  }
}// namespace RAT