// DQTellieProc.cc
// Contact person: Mark Stringer<ms711@sussex.ac.uk>
// See DQTellieProc.hh for more 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 <sstream>

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

using namespace std;
using namespace json;


namespace RAT {


    DQTellieProc::DQTellieProc() : DataQualityProc( "dqtellieproc")
    {
        fDQChecks = NULL;

        fTrigCheckThresh = 0;
        fPulseDelayThresh = 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;
        fPulseDelayMin = fDQChecks->GetD( "pulse_delay_min" );
        fCriteria["pulse_delay_min"] = fPulseDelayMin;
        fPulseDelayMax = fDQChecks->GetD( "pulse_delay_max" );
        fCriteria["pulse_delay_max"] = fPulseDelayMax;
        fTellieChipDelay  = fDQChecks->GetD( "tellie_chip_delay" );
        fCriteria["tellie_chip_delay"] = fTellieChipDelay;
        fPulseDelayThresh = fDQChecks->GetD( "pulse_delay_thresh" );
        fCriteria["pulse_delay_thresh"] = fPulseDelayThresh;
        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;
        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;
        fTACReflectedOffsetLower = fDQChecks->GetD("tac_reflected_peak_offset_low");
        fCriteria["tac_reflected_peak_offset_low"] = fTACReflectedOffsetLower;
        fTACReflectedOffsetUpper = fDQChecks->GetD("tac_reflected_peak_offset_high");
        fCriteria["tac_reflected_peak_offset_high"] = fTACReflectedOffsetUpper;
        fTriggerBit = fDQChecks->GetI("trigger_bit");
        fCriteria["trigger_bit_checked"] = fTriggerBit;


        //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);
            fSubRunNumbers.push_back(singleSubrunInfo.getMember("sub_run_number").getInteger());
            fFibres.push_back(singleSubrunInfo.getMember("fibre").getString());
            fTellieEvents.push_back(singleSubrunInfo.getMember("number_of_shots").getInteger());
            fTelliePulseRate.push_back(1000.0/singleSubrunInfo.getMember("pulse_separation").getReal());
        }

        //Intialize the vectors
        fTellieCount = std::vector<int>(fSubRunNumbers.size(),0);
        fFailPulseDelayCount = std::vector<int> (fSubRunNumbers.size(),0);
        fAvgNHit = std::vector<double> (fSubRunNumbers.size(),0.0);

        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);
        fTH2DPsupProjection  = 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);
        fPulseDelayEfficiency = 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);
        fPulseDelayEfficiencyCheck = 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());
        fPrevEventTime = 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(),
                    410,0,4100);
            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(),
                    410,0,4100);
            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(),
                    410,0,4100);
            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(),
                    1000,-500, 500);
            fTH1DCalibTime[i]->SetDirectory(0);
            //For post run monitoring
            title = " PSUP Hit Map";
            titleString = "psup_projection"+namestream.str();
            fTH2DPsupProjection[i] = new TH2D(titleString.c_str(),
                    title.c_str(),
                    300, 0, 1, 300, 0, 1);
            fTH2DPsupProjection[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 the fibre which was firing during this subrun
        int subRunNumber = ent.GetSubRunID();
        int subRunIndex = -1000;
        for(unsigned int i=0; i<fSubRunNumbers.size(); i++){
            if(fSubRunNumbers[i] == subRunNumber){
                subRunIndex = i;
                break;
            }
        }

        if(subRunIndex == -1000){
            detail << "Event is not in subrun range from couch DB subrun number is: "<<subRunNumber<<newline;
            detail << "Sub runs in couchDB: [";
            for(unsigned int i=0; i<fSubRunNumbers.size(); i++){
                detail << fSubRunNumbers[i]<<", ";
            }
            detail << "]"<<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) )
        {
            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();
        }

        //TELLIE triggered checks
        if ( passedTriggerCheck == true)
        {
            //Check 3 - Check pulse delay
            if( fFirstTellieEventTimeSet[subRunIndex] == false && ev.GetGTID() != 0 )
            {
                fPrevEventTime[subRunIndex] = ev.GetUniversalTime();
                fFirstTellieEventTimeSet[subRunIndex] = true;
            }
            else if(ev.GetGTID() != 0){
                fCurEventTime[subRunIndex] = fLastEventTime[subRunIndex];
                DS::UniversalTime utDelay  = fCurEventTime[subRunIndex] - fPrevEventTime[subRunIndex];
                double pulseDelay = utDelay.GetSeconds()+(utDelay.GetNanoSeconds()*1e-9);
                fPrevEventTime = fCurEventTime;
                double lowerPulseSeparation = 1.0/fTelliePulseRate[subRunIndex] + fPulseDelayMin + fTellieChipDelay; // delay + upper limit + chip delay (const)
                double upperPulseSeparation = 1.0/fTelliePulseRate[subRunIndex] + fPulseDelayMax + fTellieChipDelay;


                if ( pulseDelay < lowerPulseSeparation || pulseDelay > upperPulseSeparation )
                {
                    fFailPulseDelayCount[subRunIndex]++;
                }
            }
            const int NHit = ev.GetNhits();
            fTH1INHits[subRunIndex]->Fill(NHit);

            //Check 4 - NHit
            if ( NHit > fMaxNHit )
            {
                fFailMaxNHitCount[subRunIndex]++;
            }

            //Check 5 - Average 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 << " is at 0,0,0! Skipping" << newline;
                    continue;
                }
                TVector2 projection = IcosProject(pmtPos);
                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;
    }

    bool DQTellieProc::CheckNeighbours(int subRunIndex , int maxPMT)
    {
        const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
        bool isHot = false;
        int neighbour1 = maxPMT+1;
        int neighbour2 = maxPMT+2;
        int pmtCount = pmtInfo.GetCount();
        //Check to see if pmt number is valid
        if (neighbour1 >= pmtCount)
        {
            neighbour1 = 0;
            neighbour2 = 1;
        }
        else if (neighbour2 >= pmtCount)
        {
            neighbour2 = 0;
        }
        //Check neighbour1 is on same panel as maxPMT
        while (pmtInfo.GetPanelNumber(maxPMT) != pmtInfo.GetPanelNumber(neighbour1))
        {
            neighbour1 += 1;
            if (neighbour1 >= pmtCount)
            {
                neighbour1 = 0;
            }
        }
        //Check neighbour2 is on same panel as maxPMT and not neighbour1
        while ((pmtInfo.GetPanelNumber(maxPMT) != pmtInfo.GetPanelNumber(neighbour2)) ||  (neighbour1 == neighbour2))
        {
            neighbour2 += 1;
            if (neighbour2 >= pmtCount)
            {
                neighbour2 = 0;
            }
        }
        // If both neighbour 1 and neighbour 2 has a hit count < 80% of maxPMT then it is assumed that maxPMT is hot
        if ( fHitCount[subRunIndex][neighbour1] < 0.8*fHitCount[subRunIndex][maxPMT] && fHitCount[subRunIndex][neighbour2] < 0.8*fHitCount[subRunIndex][maxPMT])
        {
            isHot = true;
            fTH1DLcnHits[subRunIndex]->SetBinContent(maxPMT,0);
        }
        return isHot;
    }

    void DQTellieProc::PlotFibreCheck(const TVector2& directProject,const TVector2& reflectedProject, const TVector2& lightOriginProject, int subRunIndex)
    {
        TCanvas *C1 = new TCanvas("C1");
        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->SetLogz();
        fTH2DPsupProjection[subRunIndex]->Draw("colz");
        TMarker directSpot;
        TMarker reflectedSpot;
        TMarker lightOriginSpot;
        directSpot.SetMarkerStyle(2);
        reflectedSpot.SetMarkerStyle(3);
        lightOriginSpot.SetMarkerStyle(4);
        TLegend * leg = new TLegend(0,0,1,1);
        leg->AddEntry(&directSpot,"Direct Light Spot","p");
        leg->AddEntry(&reflectedSpot,"Reflected Light Spot","p");
        leg->AddEntry(&lightOriginSpot,"Calculated Light Origin","p");
        directSpot.DrawMarker(directProject.X(),directProject.Y());
        reflectedSpot.DrawMarker(reflectedProject.X(),reflectedProject.Y());
        lightOriginSpot.DrawMarker(lightOriginProject.X(),lightOriginProject.Y());
        legendPad->cd();
        leg->Draw();
        C1->cd();
        legendPad->Draw();
        histoPad->Draw();
        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 leg;
        leg = NULL;
        delete histoPad;
        histoPad = NULL;
        delete legendPad;
        legendPad = NULL;
        delete C1;
        C1 = NULL;
    }

    void DQTellieProc::FibreCheck(int subRunIndex)
    {
        // Get the most frequently hit PMT
        bool isHot = true;
        int maxPMT = -9999;
        const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
        WriteToRoot(fTH1DLcnHits[subRunIndex]); //Write before manip
        while(isHot == true)
        {
            maxPMT = fTH1DLcnHits[subRunIndex]->GetMaximumBin();
            isHot = CheckNeighbours(subRunIndex, maxPMT);
        }
        TVector3 maxPMTPos = pmtInfo.GetPosition(maxPMT);
        int numDirect = 0;
        int numReflected = 0;
        double posDirect[3] = {0.0, 0.0, 0.0};
        double posReflected[3] = {0.0, 0.0, 0.0};

        debug << "DQTellieProc::FibreCheck: Max PMT Pos x: " << maxPMTPos.X() << " y " << maxPMTPos.Y() << " z " << maxPMTPos.Z() << endl;

        for(unsigned int i=0;i<fHitCount[subRunIndex].size();i++)
        {
            if(pmtInfo.GetType(i) != DU::PMTInfo::NORMAL || fHitCount[subRunIndex][i]==0)
                continue;

            TVector3 pmtPos = pmtInfo.GetPosition(i);
            if(pmtPos.Dot(maxPMTPos)> 0)
            {
                // direct
                numDirect += fHitCount[subRunIndex][i];
                posDirect[0] += fHitCount[subRunIndex][i] * pmtPos.X();
                posDirect[1] += fHitCount[subRunIndex][i] * pmtPos.Y();
                posDirect[2] += fHitCount[subRunIndex][i] * pmtPos.Z();
            }
            else if(pmtPos.Dot(maxPMTPos) < 0)
            {
                //reflected
                numReflected += fHitCount[subRunIndex][i];
                posReflected[0] += fHitCount[subRunIndex][i] * pmtPos.X();
                posReflected[1] += fHitCount[subRunIndex][i] * pmtPos.Y();
                posReflected[2] += fHitCount[subRunIndex][i] * pmtPos.Z();
            }
        }
        double radDirect = 0;
        double radReflected = 0;
        for(int i=0;i<3;i++)
        {
            posDirect[i] /= static_cast<double>(numDirect);
            posReflected[i] /= static_cast<double>(numReflected);
            radDirect += (posDirect[i] * posDirect[i]);
            radReflected += (posReflected[i] * posReflected[i]);
        }
        radDirect = sqrt(radDirect);
        radReflected = sqrt(radReflected);

        TVector3 directLight;
        TVector3 reflectedLight;

        directLight.SetX(posDirect[0]);
        directLight.SetY(posDirect[1]);
        directLight.SetZ(posDirect[2]);

        reflectedLight.SetX(posReflected[0]);
        reflectedLight.SetY(posReflected[1]);
        reflectedLight.SetZ(posReflected[2]);

        TVector3 lightOrigin;
        double psupRadius = 8406.75;
        lightOrigin = psupRadius*(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());

        detail << "DQTellieProc::FibreCheck: Direct light centre: x: " << posDirect[0] << " y: " << posDirect[1] << " z: " << posDirect[2] << newline;
        detail << "Reflected light centre: x: " << posReflected[0] << " y: " << posReflected[1] << " z: " << posReflected[2] << newline;
        detail  << "Light Origin: x: " << lightOrigin.X() << " y: " << lightOrigin.Y() << " z: " << lightOrigin.Z() << newline;

        double minDist = 100000;
        string likelyFibre;
        // There are 92 distinct fibre positions on PSUP. Fibres > 92 are spares and are in the same position as other fibres.
        int maxNumFibre = 93;
        for ( int i = 0; i < maxNumFibre; i++)
        {
            stringstream fibNum;
            fibNum.str("");
            //zero-padding fibre number
            fibNum << "FT" << setw(3) << setfill('0') << i << "A";
            DBLinkPtr fibre;
            TVector3 fibrePos;
            try{
                //Get fibre from DB
                fibre = DB::Get()->GetLink("FIBRE", fibNum.str());
                fibrePos.SetX(fibre->GetD("x"));
                fibrePos.SetY(fibre->GetD("y"));
                fibrePos.SetZ(fibre->GetD("z"));
            }
            catch ( DBNotFoundError &e ){
                warn << "WARNING: DQTellieProc::FibreCheck " << fibNum.str() << " does not exist. Skipping." << newline;
                continue;
            }
            debug << "Fibre Pos: x: " << fibrePos.X() << " y: " << fibrePos.Y() << " z: " << fibrePos.Z() << newline;
            double dist = (lightOrigin-fibrePos).Mag();
            debug << "Dist from  fibre: "<<dist<<newline;
            if (dist < minDist)
            {
                minDist = dist;
                likelyFibre = fibNum.str();
                debug << "Likely fibre: "<<likelyFibre<<newline;
            }
        }
        TVector2 directProject = IcosProject(directLight);
        TVector2 reflectedProject = IcosProject(reflectedLight);
        TVector2 lightOriginProject = IcosProject(lightOrigin);
        PlotFibreCheck(directProject,reflectedProject,lightOriginProject, subRunIndex);
        debug << "DQTellieProc::FibreCheck: Likely fibre: " << likelyFibre << " distance from calculated origin " << minDist << newline;
        if (likelyFibre == fFibres[subRunIndex])
        {
            detail << "DQTellieProc::FibreCheck: run " << GetRunID() << "subrun "<<fSubRunNumbers[subRunIndex]<<" passed the fibre check (the correct fibre is firing)." << newline;
            fCheckResult = true;
        }
        else
        {
            detail << "DQTellieProc::FibreCheck: run " << GetRunID() << "subrun "<<fSubRunNumbers[subRunIndex]<< " failed the fibre check." << newline;
            detail << "Fibre that should be firing: " << fFibres[subRunIndex] << ". Fibre that appears to be firing: " << likelyFibre << newline;
            fCheckResult = false;
        }
        fFibreFiringGuess[subRunIndex] = likelyFibre;
        WriteToRoot(fTH2DPsupProjection[subRunIndex]);
        if(fCheckResult){
            fCorrectFibreCheck[subRunIndex] = 1;
        }
        else{
            fCorrectFibreCheck[subRunIndex] = 0;
        }
    }

    void DQTellieProc::PlotTac(std::vector<float> xPeaks,int subRunIndex)
    {
        TCanvas *C1 = new TCanvas("C1");
        C1->cd();
        fTH1DTacCount[subRunIndex]->Draw();
        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);
            TMarker marker;
            marker.SetMarkerStyle(3);
            marker.DrawMarker(xp,yp);
        }
        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::MinMaxIndex(const std::vector<float> v, size_t& min_idx, size_t& max_idx)
    {
        if(v.size()==0){
            return;
        }
        float min = v[0];
        float max = v[0];
        size_t idx = 0;
        for(vector<float>::const_iterator it = v.begin(); it != v.end(); it++)
        {
            if(*it < min)
            {
                min = *it;
                min_idx = idx;
            }
            if(*it > max)
            {
                max = *it;
                max_idx = idx;
            }
            idx += 1;
        }
    }
    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
        TSpectrum *peakFinderRef = new TSpectrum(2*fNPeaks[subRunIndex]);
        double threshold = 0.03;
        int nFoundRef = peakFinderRef->Search(fTH1DTacCountReflected[subRunIndex],
                7., "new", threshold);
        TSpectrum *peakFinderDir = new TSpectrum(2*fNPeaks[subRunIndex]);
        threshold = 0.2;
        int nFoundDir = peakFinderDir->Search(fTH1DTacCountDirect[subRunIndex],
                7., "new", threshold);
        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;
        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);
            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);
            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 late_idx = -1;
        size_t pre_idx = -1;
        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() <<" subrun "<< fSubRunNumbers[subRunIndex]  <<" passed the number of timing peaks check (3 peaks expected and found)." << newline;
            fCheckResult = true;
        }
        else
        {
            detail << "DQTellieProc::TimingChecks: run " << GetRunID() << " subrun "<< fSubRunNumbers[subRunIndex] <<" failed the number of timing peaks check (3 peaks expected and " << fNPeaks[subRunIndex] << " 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() << "subrun " << fSubRunNumbers[subRunIndex]<< " passed the timing peak amplitude check (prompt peak at least "<<fTACPromptHeightRatio<<" times bigger)." << newline;
            fCheckResult = true;
        }
        else
        {
            detail << "DQTellieProc::TimingChecks: run " << GetRunID() << "subrun " << fSubRunNumbers[subRunIndex]<< " failed the timing peak amplitude check (Prompt peak at least "<<fTACPromptHeightRatio<<" times bigger)." << newline;
            detail << " Prompt amplitude: " << amp_prompt[0] << " pre pulse amplitude " << amp_pre[0] << " 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() << " subrun "<< fSubRunNumbers[subRunIndex]  <<" passed the prompt peak time check (between "<<fTACPromptPeakLowerLimit<<" and "<<fTACPromptPeakUpperLimit<<" ADC)." << newline;
            fCheckResult = true;
        }
        else{

            detail << "DQTellieProc::TimingChecks: run " << GetRunID() << " subrun "<< fSubRunNumbers[subRunIndex]  <<" failed the prompt peak time check (between "<<fTACPromptPeakLowerLimit<<" and "<<fTACPromptPeakUpperLimit<<" ADC)." << newline;
            detail << "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]-fTACReflectedOffsetUpper) && tac_late < (tac_prompt[0]-fTACReflectedOffsetLower) && tac_pre > (tac_prompt[0]+fTACReflectedOffsetLower) && tac_pre[0] < (tac_prompt[0]+fTACReflectedOffsetUpper))
        {
            info << "DQTellieProc::TimingChecks: run " << GetRunID() << " subrun "<< fSubRunNumbers[subRunIndex]  <<" passed the peak time check. Pre pulse between "<<fTACReflectedOffsetLower<<" and "<<fTACReflectedOffsetUpper<<" ADC more than prompt. Late pulse between "<<fTACReflectedOffsetLower<<" and "<<fTACReflectedOffsetUpper<<" ADC less than prompt." << newline;
            info << " Prompt time: " << tac_prompt[0] << " ADC. Pre pulse time " << tac_pre[0] << " ADC. Late pulse time " << tac_late[0] << " ADC." << newline;
            fCheckResult = true;
        }
        else
        {
            info << "DQTellieProc::TimingChecks: run " << GetRunID() << " subrun "<< fSubRunNumbers[subRunIndex]  <<" failed the peak time check. Pre pulse between "<<fTACReflectedOffsetLower<<" and "<<fTACReflectedOffsetUpper<<" ADC more than prompt. Late pulse between "<<fTACReflectedOffsetLower<<" and "<<fTACReflectedOffsetUpper<<" ADC less than prompt." << newline;
            info << " Prompt time: " << tac_prompt[0] << " ADC. Pre pulse time " << tac_pre[0] << " ADC. Late pulse time " << tac_late[0] << " ADC." << newline;
            fCheckResult = false;
        }

        if(fCheckResult){
            fPeakTimeSpacingCheck[subRunIndex] = 1;
        }
        else{
            fPeakTimeSpacingCheck[subRunIndex] = 0;
        }
        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]);
        TCanvas *C1 = new TCanvas("C1");
        fTH1DCalibTime[subRunIndex]->SetXTitle("Time (ns)");
        fTH1DCalibTime[subRunIndex]->SetYTitle("Counts");
        fTH1DCalibTime[subRunIndex]->GetXaxis()->SetTitleOffset( 1.5 );
        string filename = fTH1DCalibTime[subRunIndex]->GetName();
        filename = filename+".png";
        C1->cd();
        fTH1DCalibTime[subRunIndex]->Draw();
        fTH1DCalibTime[subRunIndex]->Write();
        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] << "  :  " << fTellieCount[s] << " tellie events " << fTellieEvents[s] << newline;
            if( trigCheckResult >= fTrigCheckThresh && trigCheckResult <= 100.0)
            {
                detail << "DQTellieProc::EndOfRun: run " << GetRunID() <<" subrun "<<fSubRunNumbers[s]<< " passed trigger check" << newline;
                fCheckResult = true;
            }
            else if( trigCheckResult > 100.0 )
            {
                detail << "DQTellieProc::EndOfRun: run " << GetRunID() << " subrun "<<fSubRunNumbers[s]<< " failed trigger check" << newline;
                detail << " --> " << trigCheckResult << "% of events had the correct active triggers. Number of tellie events in check is wrong!" << newline;
                fCheckResult = false;
            }
            else
            {
                detail << "DQRunProc::EndOfRun: run " << GetRunID() << " subrun "<<fSubRunNumbers[s] << " failed trigger check" << newline;
                detail << " --> " << trigCheckResult << "% of events (less than " << fTrigCheckThresh;
                detail << "% threshold) had the correct active triggers" << newline;
                detail << (fTellieEvents[s] - fTellieCount[s]) << " events did not trigger." << 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() << " subrun "<<fSubRunNumbers[s]<<" time from 10 Mhz Clock (last (= " << fLastEventTime[s].GetSeconds();
            debug << ") - first (= " << fFirstEventTime[s].GetSeconds() << ")) = " << runLength << newline;
            fSubRunLengths[s] = runLength;
            double minRunLength = fTellieEvents[s]/fTelliePulseRate[s];

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

            //Check 3 - pulse delay checks
            double pulseDelayCheck = 1.0 -  static_cast<double>(fFailPulseDelayCount[s])/static_cast<double>(fTellieEvents[s]);
            fPulseDelayEfficiency[s]  = pulseDelayCheck;
            if (pulseDelayCheck < fPulseDelayThresh )
            {
                detail << "DQTellieProc::EndOfRun: run " << GetRunID() <<" subrun "<<fSubRunNumbers[s] << " (pulse delay efficiency = " << pulseDelayCheck;
                detail << ") passed  max/min pulse delay check" << newline;
                fCheckResult = true;
            }
            else
            {
                detail << "DQTellieProc::EndOfRun: run " << GetRunID() << "subrun "<<fSubRunNumbers[s]<< " failed pulse delay check" << newline;
                detail << " --> pulse delay efficiency (" << pulseDelayCheck;
                detail << ") was less than minimum efficiency for min/max pulse delay (" << fPulseDelayThresh << ")" << newline;
                fCheckResult = false;
            }
            if(fCheckResult){
                fPulseDelayEfficiencyCheck[s] = 1;
            }
            else{
                fPulseDelayEfficiencyCheck[s] = 0;
            }

            //Check 4 - NHit
            if ( fFailMaxNHitCount[s] == 0)
            {
                detail << "DQTellieProc::EndOfRun: run " << GetRunID() << "subrun "<<fSubRunNumbers[s];
                detail << " passed the maximum NHit check" << newline;
                fCheckResult = true;
            }
            else
            {
                detail << "DQTellieProc::EndOfRun: run" << GetRunID() << "subrun "<<fSubRunNumbers[s] << " failed NHit check";
                detail << " --> " << fFailMaxNHitCount[s] << " events had more than the maximum  NHit " << 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() << " subrun "<<fSubRunNumbers[s] <<" (average NHit = " << fAvgNHit[s];
                detail << ") passed the maximum average NHit check" << newline;
                fCheckResult = true;
            }
            else
            {
                detail << "DQTellieProc::EndOfRun: run" << GetRunID() << " subrun "<< fSubRunNumbers[s] << " failed average NHit check";
                detail << " --> Average NHit (" << fAvgNHit[s] << ") was more than the maximum average NHit." << 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);
            //Draw out the NHit Histos
            WriteToRoot(fTH1INHits[s]);
            //Draw nhit plots as pngs
            TCanvas *C1 = new TCanvas("C1");
            C1->cd();
            stringstream namestream;
            namestream << "NHit Distribution "<< fSubRunNumbers[s];
            C1->SetName(namestream.str().c_str());
            fTH1INHits[s]->Draw();
            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_delay_efficiency",fPulseDelayEfficiencyCheck);
        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("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_delay_efficiency_check",fPulseDelayEfficiencyCheck);
        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 pulseDelayEfficiencyRunCheck = 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(pulseDelayEfficiencyRunCheck){
                if(fPulseDelayEfficiencyCheck[s] == 0){
                    pulseDelayEfficiencyRunCheck = 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_delay", pulseDelayEfficiencyRunCheck);
        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 fTH2DPsupProjection[s];
            fTH2DPsupProjection[s] = NULL;
            delete fTH1INHits[s];
            fTH1INHits[s] = NULL;
        }
        DataQualityProc::EndOfRun( run );
    }
}// namespace RAT