#include <RAT/DQSmellieProc.hh>
#include <RAT/DB.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/Log.hh>
using namespace RAT;
using namespace RAT::DS;

#include <string>
#include <iostream>
#include <vector>
#include <sstream>
#include <fstream>
#include <cstdlib>
using namespace std;

#include <TVector3.h>
#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TFile.h>
#include <TSpectrum.h>
#include <TMarker.h>
using namespace ROOT;
//For pi
#include <CLHEP/Units/PhysicalConstants.h>
using namespace CLHEP;


DQSmellieProc::DQSmellieProc(): DataQualityProc("DQSmellieProc")
{
}


void DQSmellieProc::BeginOfRun(DS::Run& run)
{
    DataQualityProc::BeginOfRun(run);
    fRunNumber = run.GetRunID();
    fPMTInformation = DU::Utility::Get()->GetPMTInfo();
    // Get the processor-specific information from the "General_Info" entry in the SMELLIE_MONITOR_PROC.ratdb file
    DBLinkPtr dbLink = DB::Get()->GetLink("DQCHECKS", "smellie");
    //Get run information
    DBLinkPtr runInfo = DB::Get()->GetLink("SMELLIE_RUN","");
    DBLinkPtr subRunInfo;

    fNumberOfNodes = dbLink->GetI("numberOfNodes");
    fFibreIDs = dbLink->GetSArray("fibres");
    fFirstTimePeak_Low = dbLink->GetD("firstTimePeak_LowBound");
    fFirstTimePeak_High = dbLink->GetD("firstTimePeak_HighBound");
    fSecondTimePeak_Low = dbLink->GetD("secondTimePeak_LowBound");
    fSecondTimePeak_High = dbLink->GetD("secondTimePeak_HighBound");
    fSMELLIEDigitiserChannel = dbLink->GetI("digitiserChannel");
    fDigitiserSampleNumber = dbLink->GetI("digitiserSampleNumber");
    fErrorBound = dbLink->GetD("error_bound");
    fSuperKWavelengths = dbLink->GetDArray("superKWavelength");
    fSuperKNHits = dbLink->GetDArray("superKNHits");
    fRunTriggerCutBit = dbLink->GetI("trigger_bit");






    // Get the fibre information from the SMELLIE.ratdb file
    for (unsigned int f = 0; f < fFibreIDs.size(); f++)
    {
        DBLinkPtr dbNAVLink = DB::Get()->GetLink("FIBRE", fFibreIDs[f]);

        double fibreXPosition, fibreYPosition, fibreZPosition, fibreXDirection, fibreYDirection, fibreZDirection;
        try
        {
            fibreXPosition = dbNAVLink->GetD("x");
            fibreYPosition = dbNAVLink->GetD("y");
            fibreZPosition = dbNAVLink->GetD("z");
            fibreXDirection = dbNAVLink->GetD("u");
            fibreYDirection = dbNAVLink->GetD("v");
            fibreZDirection = dbNAVLink->GetD("w");
        }
        catch (RAT::DBNotFoundError &e) {continue;}

        // Use the positions of the 0-degree fibres as the positions of the respective nodes
        if ((fFibreIDs[f] == "FS007") || (fFibreIDs[f] == "FS025") || (fFibreIDs[f] == "FS137") || (fFibreIDs[f] == "FS255"))
        {
            TVector3 nodePosition(fibreXPosition, fibreYPosition, fibreZPosition);
            fNodePositions.push_back(nodePosition);
        }

        TVector3 fibreDirection(fibreXDirection, fibreYDirection, fibreZDirection);
        fFibreDirections.push_back(fibreDirection);
    }

    // Calculate the PMT Radius as the mean radius of the 4 nodes
    double sumOfRadii = 0.0;
    for (int n = 0; n < fNumberOfNodes; n++)
    {
        double nodeRadius = sqrt((fNodePositions[n].X() * fNodePositions[n].X()) + (fNodePositions[n].Y() * fNodePositions[n].Y()) + (fNodePositions[n].Z() * fNodePositions[n].Z()));
        sumOfRadii += nodeRadius;
    }
    fPMTRadius = (sumOfRadii / fNumberOfNodes);

    // Get the run-specific information from the ratdb configuration file, and set the trigger type cut value based on the run type
    //Get smellie info from the ratdb file
    stringstream runString;
    runString << fRunNumber;


    //Get JSON array from run information file
    json::Value subrunInfo = runInfo->GetJSON("sub_run_info");
    fNumberOfSubruns = subrunInfo.getArraySize();;
    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());
        fSubrunFibreIDsVector.push_back(singleSubrunInfo.getMember("fibre").getString());
        fSubrunNumberOfEventsVector.push_back(singleSubrunInfo.getMember("number_of_shots").getInteger());
        fLaserType.push_back(singleSubrunInfo.getMember("laser").getString());
        fSubrunLaserFrequenciesVector.push_back(singleSubrunInfo.getMember("pulse_rate").getInteger());
        //If running the superK laser calculate the mean
        if(fLaserType.at(i) == "superK"){
            double lowEdge = singleSubrunInfo.getMember("wavelength_low_edge").getInteger();
            double highEdge = singleSubrunInfo.getMember("wavelength_high_edge").getInteger();
            fSubrunWavelengthsVector.push_back((lowEdge+highEdge)/2.0);
        }
        //Obtain the wavelengths from the laser name in the SMELLE_RUN file fixed wavelength lasers have the name PQ440 for 440 fixed wavelength
        else{
            double laserWl =  atof(fLaserType.at(i).substr(2,fLaserType.at(i).size()-2).c_str());
            fSubrunWavelengthsVector.push_back(laserWl);
            if(laserWl == 0.0){
                Log::Die("SMELLIEDQProc: Unable to parse fixed laser wavelength from table name");
            }

        }
        fSubrunLaserIntensitiesVector.push_back(singleSubrunInfo.getMember("intensity").getInteger());
    }

    //Setting all subrun checks to True
    smellieFibreCheck = true;
    smellieFrequencyCheck = true;
    smellieIntensityCheck = true;
    smellieNumberOfEvents = true;


}

void DQSmellieProc::EndOfRun(DS::Run& run)
{
    vector<string> subrunEntriesVector;
    // Run the checks on all subruns simultaneously
    // The check (number 6) on CAEN waveforms only outputs plots, so will be disabled if the user has set the processor to NOT draw plots
    vector< vector<double> > allSubruns_nhitsInformation = Check_1_HitInformation();
    info << "DQSmellieProc::EndOfRun: Check_1_Fundamentals - complete" << newline;
    vector<int> allSubruns_numberOfEvents = Check_2_EventTriggering(allSubruns_nhitsInformation);
    info<< "DQSmellieProc::EndOfRun: Check_2_EventTriggering - complete" << newline;
    vector<double> allSubruns_frequencies = Check_3_LaserFrequency();
    info << "DQSmellieProc::EndOfRun: Check_3_LaserFrequency - complete" << newline;
    vector<string> allSubruns_fibreIDs = Check_4_FibreID();
    info << "DQSmellieProc::EndOfRun: Check_4_FibreID - complete" << newline;
    vector<double> allSubruns_intensities = Check_5_LaserIntensity(allSubruns_nhitsInformation, allSubruns_fibreIDs);
    info << "DQSmellieProc::EndOfRun: Check_5_LaserIntensity - complete" << newline;
    Check_6_DigitiserWaveforms();
    info << "DQSmellieProc::EndOfRun: Check_6_DigitiserWaveforms - complete" << newline;
    // For each subrun ...
    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        // ... compare the measured values of the checks with those from the configuration file
        // For the frequency and intensity, the measured and config values will not be EXACTLY the same, so consider the measured value good if between +/- [fErrorBound]% of the config value
        bool numberOfEventsCheck = false;
        if (allSubruns_numberOfEvents[s] == fSubrunNumberOfEventsVector[s]) {
            numberOfEventsCheck = true;
            info << "Passed Number of events check for subrun: "<<s<<" number of events was: "<<allSubruns_numberOfEvents[s]<<" expected: "<<fSubrunNumberOfEventsVector[s]<<newline;
            smellieNumberOfEventsSubrunVector.push_back(1);
        }
        else{
            info << "Failed Number of events check for subrun: "<<s<<" number of events was: "<<allSubruns_numberOfEvents[s]<<" expected: "<<fSubrunNumberOfEventsVector[s]<<newline;
            smellieNumberOfEventsSubrunVector.push_back(0);
        }
        if(smellieNumberOfEvents){
            smellieNumberOfEvents = numberOfEventsCheck;
        }
        bool frequencyCheck = false;
        double lowBound_ActualFrequency = (1 - (fErrorBound / 100)) * fSubrunLaserFrequenciesVector[s];
        double highBound_ActualFrequency = (1 + (fErrorBound / 100)) * fSubrunLaserFrequenciesVector[s];
        if ((allSubruns_frequencies[s] >  lowBound_ActualFrequency) && (allSubruns_frequencies[s] < highBound_ActualFrequency)) {
            frequencyCheck = true;
            info << "Passed Frequency Check for subrun: "<<s<<", frequency is: "<<allSubruns_frequencies[s]<<" limits are: ["<<lowBound_ActualFrequency<<","<<highBound_ActualFrequency<<"]"<<newline;
            smellieFrequencyCheckSubrunVector.push_back(1);
        }
        else{
            info << "Failed  Frequency Check for subrun: "<<s<<", frequency is: "<<allSubruns_frequencies[s]<<" limits are: ["<<lowBound_ActualFrequency<<","<<highBound_ActualFrequency<<"]"<<newline;
            smellieFrequencyCheckSubrunVector.push_back(0);
        }
        if(smellieFrequencyCheck){
            smellieFrequencyCheck = frequencyCheck;
        }
        bool fibreIDCheck = false;
        if (allSubruns_fibreIDs[s] == fSubrunFibreIDsVector[s]) {
            fibreIDCheck = true;
            info << "Passed Fibre Check for subrun: "<<s<<", fibre firing is: "<<allSubruns_fibreIDs[s]<<" expected fibre: "<<fSubrunFibreIDsVector[s]<<newline;
            smellieFibreCheckSubrunVector.push_back(1);
        }
        else{
            info << "Failed Fibre Check for subrun: "<<s<<", fibre firing is: "<<allSubruns_fibreIDs[s]<<" expected fibre: "<<fSubrunFibreIDsVector[s]<<newline;
            smellieFibreCheckSubrunVector.push_back(0);
        }
        if(smellieFibreCheck){
            smellieFibreCheck = fibreIDCheck;
        }
        bool intensityCheck = false;
        //If using a fixed wavelength laser compare intensities
        if(fLaserType[s] != "superK"){
            double lowBound_ActualIntensity = (1 - (fErrorBound / 100)) * fSubrunLaserIntensitiesVector[s];
            double highBound_ActualIntensity = (1 + (fErrorBound / 100)) * fSubrunLaserIntensitiesVector[s];
            if ((allSubruns_intensities[s] >  lowBound_ActualIntensity) && (allSubruns_intensities[s] < highBound_ActualIntensity)) {
                intensityCheck = true;
                info << "Passed Intensity Check for subrun: "<<s<<", intensity is: "<<allSubruns_intensities[s]<<" limits are: ["<<lowBound_ActualIntensity<<","<<highBound_ActualIntensity<<"]"<<newline;
                smellieIntensityCheckSubrunVector.push_back(1);
            }
            else{
                info << "Failed Intensity Check for subrun: "<<s<<", intensity is: "<<allSubruns_intensities[s]<<" limits are: ["<<lowBound_ActualIntensity<<","<<highBound_ActualIntensity<<"]"<<newline;
                smellieIntensityCheckSubrunVector.push_back(0);
            }
        }
        else{
            double lowBound_ActualWavelength = (1 - (fErrorBound / 100)) * fSubrunWavelengthsVector[s];
            double highBound_ActualWavelength = (1 + (fErrorBound / 100)) * fSubrunWavelengthsVector[s];
            //When the Super K laser is being used the vector allSubruns_intensities  stores the guess for the lasers wavelength
            if ((allSubruns_intensities[s] >  lowBound_ActualWavelength) && (allSubruns_intensities[s] < highBound_ActualWavelength)) {
                intensityCheck = true;
                info << "Passed Intensity Check (Using SuperK laser) for subrun: "<<s<<", wavelength is: "<<allSubruns_intensities[s]<<" limits are: ["<<lowBound_ActualWavelength<<","<<highBound_ActualWavelength<<"]"<<newline;
                smellieIntensityCheckSubrunVector.push_back(1);
            }
            else{
                info << "Failed  Intensity Check (Using SuperK laser) for subrun: "<<s<<", wavelength is: "<<allSubruns_intensities[s]<<" limits are: ["<<lowBound_ActualWavelength<<","<<highBound_ActualWavelength<<"]"<<newline;
                smellieIntensityCheckSubrunVector.push_back(0);
            }

        }
        if(smellieIntensityCheck){
            smellieIntensityCheck = intensityCheck;
        }
    }


    //Setting pass fail variables in the database
    Update("smellieIntensityCheck",smellieIntensityCheck);
    Update("smellieFrequencyCheck",smellieFrequencyCheck);
    Update("smellieCorrectFibre",smellieFibreCheck);
    Update("smellieNumberOfEventsCheck",smellieNumberOfEvents);

    //Adding vectors of check results for subrun to database
    AddToRATDB("smellieIntensityCheckSubrun",smellieIntensityCheckSubrunVector);
    AddToRATDB("smellieFrequencyCheckSubrun",smellieFrequencyCheckSubrunVector);
    AddToRATDB("smellieFibreCheckSubrun",smellieFibreCheckSubrunVector);
    AddToRATDB("smellieNumberOfEventsCheckSubrun",smellieNumberOfEventsSubrunVector);

    AddToRATDB("events_passing_nhit_and_trigger",passedNHitPassedTrigger);
    AddToRATDB("events_passing_nhit_failing_trigger",passedNHitFailedTrigger);
    AddToRATDB("events_failing_nhit_passing_trigger",failedNHitPassedTrigger);
    AddToRATDB("events_failing_nhit_failing_trigger",failedNHitFailedTrigger);

    AddToRATDB("nhit_event_next_to_trigger_event",NHitEventAdjacentToTrigger);
    AddToRATDB("nhit_event_no_adjacent_trigger",NHitEventsMissingTrigger);
    AddToRATDB("trigger_event_no_adjacent_nhit",TriggerEventsMissingNhits);

    AddToRATDB("subrun_numbers",fSubRunNumbers);

    AddToRATDB("mean_nhit_smellie_trigger_subrun",subRunMeanNhit);

    AddToRATDB("number_events_expected_subrun",fSubrunNumberOfEventsVector);

    AddToRATDB("frequency_expected_subrun",fSubrunLaserFrequenciesVector);
    AddToRATDB("frequency_actual_subrun",allSubruns_frequencies);

    AddToRATDB("fibre_expected_subrun",fSubrunFibreIDsVector);
    AddToRATDB("fibre_calculated_subrun",allSubruns_fibreIDs);

    AddToRATDB("laser_type",fLaserType);
    AddToRATDB("laser_wavelengths",fSubrunWavelengthsVector);
    AddToRATDB("laser_intensities",fSubrunLaserIntensitiesVector);
    AddToRATDB("fixed_intensities_and_superk_wavelengths",allSubruns_intensities);

    DataQualityProc::EndOfRun(run);
    // Clear ALL variables and vectors
    allSubruns_nhitsInformation.clear();
    allSubruns_numberOfEvents.clear();
    allSubruns_frequencies.clear();
    allSubruns_fibreIDs.clear();
    allSubruns_intensities.clear();

    fSubrunWavelengthsVector.clear();
    fSubrunNumberOfEventsVector.clear();
    fSubrunLaserFrequenciesVector.clear();
    fSubrunFibreIDsVector.clear();
    fSubrunLaserIntensitiesVector.clear();

    fAllEvents_SubRunIDs.clear();
    fAllEvents_Nhits.clear();
    fAllEvents_TrigTypes.clear();
    fAllEvents_TickCounts.clear();
    fAllEvents_Waveforms.clear();
    fAllEvents_pmtTimes.clear();
    fAllEvents_pmtThetas.clear();
    fAllEvents_pmtPhis.clear();
}


DataQualityProc::Result DQSmellieProc::DSEvent(DS::Run&, DS::Entry& ds)
{
    // Add the values from all triggered events in this DS entry into the vectors of useful information
    for (size_t eventEVID = 0; eventEVID < ds.GetEVCount(); eventEVID++)
    {
        DS::EV& triggeredEvent = ds.GetEV(eventEVID);
        int eventSubRunID = ds.GetSubRunID();
        int eventTriggerType = triggeredEvent.GetTrigType();
        int eventNhits = triggeredEvent.GetCalPMTs().GetCount();
        ULong64_t eventTickCount = triggeredEvent.GetClockCount50();

        DS::Digitiser eventDigitiser;
        try {
            eventDigitiser = triggeredEvent.GetDigitiser();
            if (eventDigitiser.ExistsWaveform(fSMELLIEDigitiserChannel))
            {
                fAllEvents_Waveforms.push_back(eventDigitiser.GetWaveform(fSMELLIEDigitiserChannel));
            }
            else
            {
                info << "DQSmellieProc::DSEvent: SubRun " << eventSubRunID << " | Event " << eventEVID << " - no digitiser waveform exists on Channel " << fSMELLIEDigitiserChannel << newline;
                vector<UShort_t> outVector;
                outVector.push_back(65535);
                fAllEvents_Waveforms.push_back(outVector);
            }
        }
        catch (DS::DataNotFound& e) {
            info << "DQSmellieProc::DSEvent: SubRun " << eventSubRunID << " | Event " << eventEVID << " - Unable to load digitiser" << newline;
            vector<UShort_t> outVector;
            outVector.push_back(65535);
            fAllEvents_Waveforms.push_back(outVector);
        }

        fAllEvents_SubRunIDs.push_back(eventSubRunID);
        fAllEvents_TrigTypes.push_back(eventTriggerType);
        fAllEvents_Nhits.push_back(eventNhits);
        fAllEvents_TickCounts.push_back(eventTickCount);

        vector<double> eventPMTTimes, eventPMTThetas, eventPMTPhis;
        for (size_t k = 0; k < triggeredEvent.GetCalPMTs().GetCount(); k++)
        {
            const DS::PMTCal& pmt = triggeredEvent.GetCalPMTs().GetPMT(k);
            TVector3 pmtPosition = fPMTInformation.GetPosition(pmt.GetID());

            eventPMTTimes.push_back(pmt.GetTime());
            eventPMTThetas.push_back(pmtPosition.Theta());
            eventPMTPhis.push_back(pmtPosition.Phi());
        }
        fAllEvents_pmtTimes.push_back(eventPMTTimes);
        fAllEvents_pmtThetas.push_back(eventPMTThetas);
        fAllEvents_pmtPhis.push_back(eventPMTPhis);
    }

    return OKTRUE;
}


vector< vector<double> > DQSmellieProc::Check_1_HitInformation()
{
    // Prepare a set of histograms for each subrun
    vector<TH1D> allSubruns_nhitsPlots_allEvents, allSubruns_nhitsPlots_trigCutPassed, allSubruns_nhitsPlots_trigCutFailed;
    vector<TH1D> allSubruns_timePlots_allEvents, allSubruns_timePlots_trigCutPassed, allSubruns_timePlots_trigCutFailed;
    //Store flipped histograms for display in plots
    vector<TH2D> allSubruns_thetaPhiPlots_allEvents_pi_minus_theta, allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta, allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        TH1D singleSubrun_nhitsPlot_allEvents ("singleSubrun_nhitsPlot_allEvents", "Nhits per Event for All Events; Nhits", 250, 0.0, 500.0);
        allSubruns_nhitsPlots_allEvents.push_back(singleSubrun_nhitsPlot_allEvents);
        TH1D singleSubrun_nhitsPlot_trigCutPassed ("singleSubrun_nhitsPlot_trigCutPassed", "Nhits per Event for Events Passing the Trigger Cut; Nhits", 250, 0.0, 500.0);
        allSubruns_nhitsPlots_trigCutPassed.push_back(singleSubrun_nhitsPlot_trigCutPassed);
        TH1D singleSubrun_nhitsPlot_trigCutFailed ("singleSubrun_nhitsPlot_trigCutFailed", "Nhits per Event for Events Failing the Trigger Cut; Nhits", 250, 0.0, 500.0);
        allSubruns_nhitsPlots_trigCutFailed.push_back(singleSubrun_nhitsPlot_trigCutFailed);

        TH1D singleSubrun_timePlot_allEvents ("singleSubrun_timePlot_allEvents", "Time of each Triggered PMT for All Events; Time, ns", 500, 0, 500);
        allSubruns_timePlots_allEvents.push_back(singleSubrun_timePlot_allEvents);
        TH1D singleSubrun_timePlot_trigCutPassed ("singleSubrun_timePlot_trigCutPassed", "Time of each Triggered PMT for Events Passing the Trigger Cut; Time, ns", 500, 0, 500.0);
        allSubruns_timePlots_trigCutPassed.push_back(singleSubrun_timePlot_trigCutPassed);
        TH1D singleSubrun_timePlot_trigCutFailed ("singleSubrun_timePlot_trigCutFailed", "Time of each Triggered PMT for Events Failing the Trigger Cut; Time, ns", 500, 0, 500.0);
        allSubruns_timePlots_trigCutFailed.push_back(singleSubrun_timePlot_trigCutFailed);

        TH2D singleSubrun_thetaPhiPlot_allEvents_pi_minus_theta("singleSubrun_thetaPhiPlot_allEvents_pi_minus_theta", "Theta and Phi of each Triggered PMT for All Events; Phi, rads; #pi - Theta, rads", 300, -3.5, 3.5, 300, 0, 3.5);
        allSubruns_thetaPhiPlots_allEvents_pi_minus_theta.push_back(singleSubrun_thetaPhiPlot_allEvents_pi_minus_theta);

        TH2D singleSubrun_thetaPhiPlot_trigCutPassed_pi_minus_theta("singleSubrun_thetaPhiPlot_trigCutPassed_pi_minus_theta", "Theta and Phi of each Triggered PMT for Events Passing the Trigger Cut; Phi, rads; #pi-Theta, rads",300, -3.5, 3.5, 300, 0, 3.5);
        allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta.push_back(singleSubrun_thetaPhiPlot_trigCutPassed_pi_minus_theta);

        TH2D singleSubrun_thetaPhiPlot_trigCutFailed_pi_minus_theta ("singleSubrun_thetaPhiPlot_trigCutFailed_pi_minus_theta", "Theta and Phi of each Triggered PMT for Events Failing the Trigger Cut; Phi, rads; #pi-Theta, rads", 300, -3.5, 3.5, 300, 0.0, 3.5);
        allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta.push_back(singleSubrun_thetaPhiPlot_trigCutFailed_pi_minus_theta);

    }
    // Loop over all of the event information, and fill the appropriate set of histograms based on the event's subrun ID and whether or not it passed the trigger cut
    for (unsigned int e = 0; e < fAllEvents_SubRunIDs.size(); e++)
    {
        int eventSubrunID = fAllEvents_SubRunIDs[e];
        if(eventSubrunID >= fNumberOfSubruns){
            info << "Subrun number: "<<eventSubrunID<<" Greater than or eqaul  Number of subruns: "<<fNumberOfSubruns<<newline;
            continue;
        }
        (allSubruns_nhitsPlots_allEvents[eventSubrunID]).Fill(fAllEvents_Nhits[e]);
        for (unsigned int p = 0; p < fAllEvents_pmtTimes[e].size(); p++)
        {
            (allSubruns_timePlots_allEvents[eventSubrunID]).Fill((fAllEvents_pmtTimes[e])[p]);
            (allSubruns_thetaPhiPlots_allEvents_pi_minus_theta[eventSubrunID]).Fill((fAllEvents_pmtPhis[e])[p],pi-(fAllEvents_pmtThetas[e])[p]);
        }
        if ((fAllEvents_TrigTypes[e] >> fRunTriggerCutBit) & 1)
        {
            (allSubruns_nhitsPlots_trigCutPassed[eventSubrunID]).Fill(fAllEvents_Nhits[e]);
            for (unsigned int p = 0; p < fAllEvents_pmtTimes[e].size(); p++)
            {

                (allSubruns_timePlots_trigCutPassed[eventSubrunID]).Fill((fAllEvents_pmtTimes[e])[p]);
                (allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta[eventSubrunID]).Fill((fAllEvents_pmtPhis[e])[p], pi-(fAllEvents_pmtThetas[e])[p]);
            }
        }
        else
        {
            (allSubruns_nhitsPlots_trigCutFailed[eventSubrunID]).Fill(fAllEvents_Nhits[e]);
            for (unsigned int p = 0; p < fAllEvents_pmtTimes[e].size(); p++)
            {
                (allSubruns_timePlots_trigCutFailed[eventSubrunID]).Fill((fAllEvents_pmtTimes[e])[p]);
                (allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta[eventSubrunID]).Fill((fAllEvents_pmtPhis[e])[p],pi-(fAllEvents_pmtThetas[e])[p]);
            }
        }
    }
    // For each subrun, and if applicable, make the plots presentable and then save to the output file
    stringstream mainPlotFileNameStream;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        (allSubruns_nhitsPlots_allEvents[s]).SetStats(0);
        (allSubruns_nhitsPlots_allEvents[s]).SetLineColor(kBlue);
        (allSubruns_nhitsPlots_allEvents[s]).SetLineWidth(2);
        (allSubruns_nhitsPlots_trigCutPassed[s]).SetStats(0);
        (allSubruns_nhitsPlots_trigCutPassed[s]).SetLineColor(kGreen+2);
        (allSubruns_nhitsPlots_trigCutPassed[s]).SetLineWidth(2);
        (allSubruns_nhitsPlots_trigCutFailed[s]).SetStats(0);
        (allSubruns_nhitsPlots_trigCutFailed[s]).SetLineColor(kRed);
        (allSubruns_nhitsPlots_trigCutFailed[s]).SetLineWidth(2);

        stringstream canvas1Namestream;
        canvas1Namestream << "SMELLIEDQ_Check1_nhitsPlotsCanvas_subrun" << s;
        TCanvas nhitsPlotsCanvas ((canvas1Namestream.str()).c_str(), "Nhits per Event", 1500, 900);
        nhitsPlotsCanvas.Divide(2, 2);
        nhitsPlotsCanvas.cd(1);
        (allSubruns_nhitsPlots_allEvents[s]).Draw();
        nhitsPlotsCanvas.cd(2);
        (allSubruns_nhitsPlots_trigCutPassed[s]).Draw();
        nhitsPlotsCanvas.cd(4);
        (allSubruns_nhitsPlots_trigCutFailed[s]).Draw();
        WriteToRoot(&nhitsPlotsCanvas);
        string filename = string(nhitsPlotsCanvas.GetName())+".png";
        nhitsPlotsCanvas.Print(filename.c_str());


        (allSubruns_timePlots_allEvents[s]).SetStats(0);
        (allSubruns_timePlots_allEvents[s]).SetLineWidth(2);
        (allSubruns_timePlots_allEvents[s]).SetLineColor(kBlue);

        (allSubruns_thetaPhiPlots_allEvents_pi_minus_theta[s]).SetStats(0);
        (allSubruns_thetaPhiPlots_allEvents_pi_minus_theta[s]).SetMarkerStyle(20);
        (allSubruns_thetaPhiPlots_allEvents_pi_minus_theta[s]).SetMarkerSize(0.1);
        (allSubruns_thetaPhiPlots_allEvents_pi_minus_theta[s]).SetMarkerColor(kBlue);

        stringstream canvas2Namestream;
        canvas2Namestream << "SMELLIEDQ_Check1_hitMaps_allEvents_subrun" << s;
        TCanvas hitMapsCanvas_allEvents ((canvas2Namestream.str()).c_str(), "Theta, Phi and Time of each Triggered PMT for All Events", 1500, 900);
        hitMapsCanvas_allEvents.Divide(2, 2);
        hitMapsCanvas_allEvents.cd(1);
        (allSubruns_timePlots_allEvents[s]).Draw();
        TVirtualPad * plotPad = hitMapsCanvas_allEvents.cd(4);
        plotPad->SetLogz();
        (allSubruns_thetaPhiPlots_allEvents_pi_minus_theta[s]).Draw("colz");
        WriteToRoot(&hitMapsCanvas_allEvents);
        filename = string(hitMapsCanvas_allEvents.GetName())+".png";
        hitMapsCanvas_allEvents.Print(filename.c_str());

        (allSubruns_timePlots_trigCutPassed[s]).SetStats(0);
        (allSubruns_timePlots_trigCutPassed[s]).SetLineWidth(2);
        (allSubruns_timePlots_trigCutPassed[s]).SetLineColor(kGreen+2);

        (allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta[s]).SetStats(0);
        (allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta[s]).SetMarkerStyle(20);
        (allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta[s]).SetMarkerSize(0.1);
        (allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta[s]).SetMarkerColor(kGreen+2);

        stringstream canvas3Namestream;
        canvas3Namestream << "SMELLIEDQ_Check1_hitMaps_trigCutPassed_subrun" << s;
        TCanvas hitMapsCanvas_trigCutPassed ((canvas3Namestream.str()).c_str(), "Theta, Phi and Time of each Triggered PMT for Events Passing the Trigger Cut", 1500, 900);
        hitMapsCanvas_trigCutPassed.Divide(2, 2);
        hitMapsCanvas_trigCutPassed.cd(1);
        (allSubruns_timePlots_trigCutPassed[s]).Draw();
        plotPad = hitMapsCanvas_trigCutPassed.cd(4);
        plotPad->SetLogz();
        (allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta[s]).Draw("colz");
        WriteToRoot(&hitMapsCanvas_trigCutPassed);
        filename = string(hitMapsCanvas_trigCutPassed.GetName())+".png";
        hitMapsCanvas_trigCutPassed.Print(filename.c_str());

        (allSubruns_timePlots_trigCutFailed[s]).SetStats(0);
        (allSubruns_timePlots_trigCutFailed[s]).SetLineWidth(2);
        (allSubruns_timePlots_trigCutFailed[s]).SetLineColor(kRed);

        (allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta[s]).SetStats(0);
        (allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta[s]).SetMarkerStyle(20);
        (allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta[s]).SetMarkerSize(0.1);
        (allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta[s]).SetMarkerColor(kRed);

        stringstream canvas4Namestream;
        canvas4Namestream << "SMELLIEDQ_Check1_hitMaps_trigCutFailed_subrun" << s;
        TCanvas hitMapsCanvas_trigCutFailed ((canvas4Namestream.str()).c_str(), "Theta, Phi and Time of each Triggered PMT for Events Failing the Trigger Cut", 1500, 900);
        hitMapsCanvas_trigCutFailed.Divide(2, 2);
        hitMapsCanvas_trigCutFailed.cd(1);
        (allSubruns_timePlots_trigCutFailed[s]).Draw();
        plotPad = hitMapsCanvas_trigCutFailed.cd(4);
        plotPad->SetLogz();
        (allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta[s]).Draw("colz");
        WriteToRoot(&hitMapsCanvas_trigCutFailed);
        filename = string(hitMapsCanvas_trigCutFailed.GetName())+".png";
        hitMapsCanvas_trigCutFailed.Print(filename.c_str());
    }


    // For each subrun, find the low-edge and the mean of the main Nhits peak for events passing the Trigger Type cut
    vector< vector<double> > allSubruns_nhitsInformation;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        double lowEdgeOfNhitsPeak = (allSubruns_nhitsPlots_trigCutPassed[s]).GetBinLowEdge((allSubruns_nhitsPlots_trigCutPassed[s]).FindFirstBinAbove(20, 1));
        lowEdgeOfNhitsPeak -= 20.0;
        vector<double> singleSubrun_nhitsInformation;
        TF1* nhitsFit = new TF1("gaus", "gaus", 0.0, 500.0);
        (allSubruns_nhitsPlots_trigCutPassed[s]).Fit(nhitsFit, "RQNO");
        double meanNhits = nhitsFit->GetParameter(1);
        delete nhitsFit;
        singleSubrun_nhitsInformation.push_back(lowEdgeOfNhitsPeak);
        singleSubrun_nhitsInformation.push_back(meanNhits);
        allSubruns_nhitsInformation.push_back(singleSubrun_nhitsInformation);
        subRunMeanNhit.push_back(meanNhits);
    }

    // Clear all vectors used in this function
    allSubruns_nhitsPlots_allEvents.clear();
    allSubruns_nhitsPlots_trigCutPassed.clear();
    allSubruns_nhitsPlots_trigCutFailed.clear();
    allSubruns_timePlots_allEvents.clear();
    allSubruns_timePlots_trigCutPassed.clear();
    allSubruns_timePlots_trigCutFailed.clear();
    allSubruns_thetaPhiPlots_allEvents_pi_minus_theta.clear();
    allSubruns_thetaPhiPlots_trigCutPassed_pi_minus_theta.clear();
    allSubruns_thetaPhiPlots_trigCutFailed_pi_minus_theta.clear();
    return allSubruns_nhitsInformation;
}


vector<int> DQSmellieProc::Check_2_EventTriggering(vector< vector<double> > allSubruns_nhitsInformation)
{
    // Prepare a histogram and vectors for each subrun
    vector<TH2D> allSubruns_nhitsVsTriggerPlots_allEvents;
    vector< vector<int> > allSubruns_eventStatus_allEvents, allSubruns_eventNhits_allEvents, allSubruns_eventTrigTypes_allEvents;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        stringstream plotNamestream;
        plotNamestream << "SMELLIEDQ_Check2_nhitsVsTriggerPlot_subrun" << s;
        TH2D singleSubrun_nhitsVsTriggerPlot ((plotNamestream.str()).c_str(), "Nhits Vs. Trigger Type for All Events; Nhits; Trigger Type", 500, 0, 500, 600, 0, 600);
        allSubruns_nhitsVsTriggerPlots_allEvents.push_back(singleSubrun_nhitsVsTriggerPlot);

        vector<int> singleSubrun_eventStatus_allEvents, singleSubrun_eventNhits_allEvents, singleSubrun_eventTrigTypes_allEvents;
        allSubruns_eventStatus_allEvents.push_back(singleSubrun_eventStatus_allEvents);
        allSubruns_eventNhits_allEvents.push_back(singleSubrun_eventNhits_allEvents);
        allSubruns_eventTrigTypes_allEvents.push_back(singleSubrun_eventTrigTypes_allEvents);
    }

    // Assign a "status" to each event, based on if they passed or failed the Nhits and Trigger Type cuts and fill the appropriate histogram and vector based on the event's subrun ID
    for (unsigned int e = 0; e < fAllEvents_SubRunIDs.size(); e++)
    {
        int eventSubrunID = fAllEvents_SubRunIDs[e];
        if(eventSubrunID >= fNumberOfSubruns){
            info << "Subrun number: "<<eventSubrunID<<" Greater than or eqaul  Number of subruns: "<<fNumberOfSubruns<<newline;
            continue;
        }

        (allSubruns_nhitsVsTriggerPlots_allEvents[eventSubrunID]).Fill(fAllEvents_Nhits[e], fAllEvents_TrigTypes[e]);

        if ((fAllEvents_Nhits[e] >= (allSubruns_nhitsInformation[eventSubrunID])[0]) && ((fAllEvents_TrigTypes[e] >> fRunTriggerCutBit) & 1)) {(allSubruns_eventStatus_allEvents[eventSubrunID]).push_back(1);}
        else if ((fAllEvents_Nhits[e] >= (allSubruns_nhitsInformation[eventSubrunID])[0]) && ((fAllEvents_TrigTypes[e] >> fRunTriggerCutBit) & 1)) {(allSubruns_eventStatus_allEvents[eventSubrunID]).push_back(2);}
        else if ((fAllEvents_Nhits[e] < (allSubruns_nhitsInformation[eventSubrunID])[0]) && ((fAllEvents_TrigTypes[e] >> fRunTriggerCutBit) & 1)) {(allSubruns_eventStatus_allEvents[eventSubrunID]).push_back(3);}
        else {(allSubruns_eventStatus_allEvents[eventSubrunID]).push_back(4);}

        (allSubruns_eventNhits_allEvents[eventSubrunID]).push_back(fAllEvents_Nhits[e]);
        (allSubruns_eventTrigTypes_allEvents[eventSubrunID]).push_back(fAllEvents_TrigTypes[e]);
    }

    // For each subrun, and if applicable, make the plot presentable and save to output file

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        (allSubruns_nhitsVsTriggerPlots_allEvents[s]).SetStats(0);
        (allSubruns_nhitsVsTriggerPlots_allEvents[s]).SetMarkerStyle(20);
        (allSubruns_nhitsVsTriggerPlots_allEvents[s]).SetMarkerSize(0.5);
        (allSubruns_nhitsVsTriggerPlots_allEvents[s]).SetMarkerColor(kBlue);

        WriteToRoot(&allSubruns_nhitsVsTriggerPlots_allEvents[s]);
        string filename = allSubruns_nhitsVsTriggerPlots_allEvents[s].GetName();
        TCanvas outCan (filename.c_str(),filename.c_str(),1500,900);
        outCan.cd();
        allSubruns_nhitsVsTriggerPlots_allEvents[s].Draw();
        filename = filename+".png";
        outCan.Print(filename.c_str());


    }

    // For each subrun, output information about how many events had each type of status and save the number of true SMELLIE events
    vector<int> allSubruns_numberOfEvents;

    stringstream textFileNameStream;
    textFileNameStream << "SMELLIEDQ" << "_OutputText_Run" << fRunNumber << ".txt";
    ofstream textOutput ((textFileNameStream.str()).c_str(), ofstream::out | ofstream::app);
    info << "Making text file" << newline;
    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        info <<"On Subrun: "<<s<<newline;
        textOutput << newline;
        textOutput << "--------------------------------------------------" << newline;
        textOutput << "Check 2: Subrun " << s << " - Event Triggering" << newline;
        textOutput << newline;

        int counter_eventStatus1 = count((allSubruns_eventStatus_allEvents[s]).begin(), (allSubruns_eventStatus_allEvents[s]).end(), 1);
        int counter_eventStatus2 = count((allSubruns_eventStatus_allEvents[s]).begin(), (allSubruns_eventStatus_allEvents[s]).end(), 2);
        int counter_eventStatus3 = count((allSubruns_eventStatus_allEvents[s]).begin(), (allSubruns_eventStatus_allEvents[s]).end(), 3);
        int counter_eventStatus4 = count((allSubruns_eventStatus_allEvents[s]).begin(), (allSubruns_eventStatus_allEvents[s]).end(), 4);

        textOutput << "A true SMELLIE event in this subrun could be identified by passing either of the following cuts:" << newline;
        textOutput << " - Nhits >= " << (allSubruns_nhitsInformation[s])[0] << "  (Nhits Cut)" << newline;
        textOutput << " - Trigger Type >= " << fRunTriggerCutBit << "  (Trigger Cut)" << newline;
        textOutput << newline;
        textOutput << "There were a total of " << (allSubruns_eventStatus_allEvents[s]).size() << " triggered events in this subrun" << newline;
        textOutput << newline;
        textOutput << "Out of the " << (allSubruns_eventStatus_allEvents[s]).size() << " total events, " << counter_eventStatus1 + counter_eventStatus2 << " passed the Nhits Cut" << newline;
        textOutput << "Out of these " << counter_eventStatus1 + counter_eventStatus2 << " events, " << counter_eventStatus1 << " passed the Trigger Cut as well" << newline;
        textOutput << "Therefore, there are " << counter_eventStatus2 << " events that have the correct Nhits for a SMELLIE event but do not have the correct Trigger Type (\"Missing Triggers\")" << newline;
        textOutput << newline;
        textOutput << "Out of the " << (allSubruns_eventStatus_allEvents[s]).size() << " total events, " << counter_eventStatus3 + counter_eventStatus4 << " failed the Nhits Cut" << newline;
        textOutput << "Out of these " << counter_eventStatus3 + counter_eventStatus4 << " events, " << counter_eventStatus4 << " failed the Trigger Cut as well" << newline;
        textOutput << "Therefore, there are " << counter_eventStatus3 << " events that do not have the correct Nhits for a SMELLIE event but have the correct Trigger Type (\"Extra Triggers\")" << newline;
        passedNHitPassedTrigger.push_back(counter_eventStatus1);
        passedNHitFailedTrigger.push_back(counter_eventStatus2);
        failedNHitPassedTrigger.push_back(counter_eventStatus3);
        failedNHitFailedTrigger.push_back(counter_eventStatus4);

        // Look more closely at those events with status 2 or 3 ...
        int counter_status2with3 = 0;
        vector<int> status2with3_previousEventStatus, status2with3_previousEventTrigger, status2with3_previousEventNhits;
        vector<int> status2with3_currentEventStatus, status2with3_currentEventTrigger, status2with3_currentEventNhits;
        vector<int> status2with3_nextEventStatus, status2with3_nextEventTrigger, status2with3_nextEventNhits;
        int counter_status2without3 = 0;
        vector<int> status2without3_currentEventStatus, status2without3_currentEventTrigger, status2without3_currentEventNhits;
        int counter_status3without2 = 0;
        vector<int> status3without2_currentEventNhits;

        for (unsigned int e = 0; e < (allSubruns_eventStatus_allEvents[s]).size(); e++)
        {
            int previousEventStatus;
            if (e == 0) {previousEventStatus = 0;}
            else {previousEventStatus = (allSubruns_eventStatus_allEvents[s])[e - 1];}

            int nextEventStatus;
            if (e == ((allSubruns_eventStatus_allEvents[s]).size() - 1)) {nextEventStatus = 0;}
            else {nextEventStatus = (allSubruns_eventStatus_allEvents[s])[e + 1];}

            // For a Status 2 event, check if it is immediately preceded or followed by a Status 3 event
            if ((allSubruns_eventStatus_allEvents[s])[e] == 2)
            {
                if  ((previousEventStatus == 3) || (nextEventStatus == 3))
                {
                    counter_status2with3++;
                    status2with3_previousEventStatus.push_back(previousEventStatus);
                    status2with3_previousEventTrigger.push_back((allSubruns_eventTrigTypes_allEvents[s])[e - 1]);
                    status2with3_previousEventNhits.push_back((allSubruns_eventNhits_allEvents[s])[e - 1]);
                    status2with3_currentEventStatus.push_back((allSubruns_eventStatus_allEvents[s])[e]);
                    status2with3_currentEventTrigger.push_back((allSubruns_eventTrigTypes_allEvents[s])[e]);
                    status2with3_currentEventNhits.push_back((allSubruns_eventNhits_allEvents[s])[e]);
                    status2with3_nextEventStatus.push_back(nextEventStatus);
                    status2with3_nextEventTrigger.push_back((allSubruns_eventTrigTypes_allEvents[s])[e + 1]);
                    status2with3_nextEventNhits.push_back((allSubruns_eventNhits_allEvents[s])[e + 1]);
                }
                else
                {
                    counter_status2without3++;
                    status2without3_currentEventStatus.push_back((allSubruns_eventStatus_allEvents[s])[e]);
                    status2without3_currentEventTrigger.push_back((allSubruns_eventTrigTypes_allEvents[s])[e]);
                    status2without3_currentEventNhits.push_back((allSubruns_eventNhits_allEvents[s])[e]);
                }
            }

            // For a Status 3 event, check if it is immediately preceded or followed by a Status 2 event
            // If NOT, save the Nhits so that a user can find the event later from the nhits vs. trigger type plot
            else if ((allSubruns_eventStatus_allEvents[s])[e] == 3)
            {
                if ((previousEventStatus == 2) || (nextEventStatus == 2)) continue;
                counter_status3without2++;
                status3without2_currentEventNhits.push_back((allSubruns_eventNhits_allEvents[s])[e]);
            }
        }
        info << "Making output info"<<newline;
        // Output information that can be used for diagnosing the Missing and Extra trigger events
        textOutput << newline;
        textOutput << "Out of the " << counter_eventStatus2 << " \"Missing Trigger\" events, " << counter_status2with3 << " events are either immediately preceded or followed by an \"Extra Trigger\" event" << newline;
        NHitEventAdjacentToTrigger.push_back(counter_status2with3);
        textOutput << "Details of these events are given below (* denotes the current event index):" << newline;
        textOutput << newline;
        textOutput << "\t" << "Status[*-1]" << "\t\t" << "Trigger[*-1]" << "\t" << "Nhits[*-1]"
            << "\t | \t" << "Status[*]" << "\t" << "Trigger[*]" << "\t" << "Nhits[*]" << "\t | \t"
            << "Status[*+1]" << "\t\t" << "Trigger[*+1]" << "\t" << "Nhits[*+1]" << newline;
        for (int b = 0; b < counter_status2with3; b++)
        {
            textOutput << "\t" << status2with3_previousEventStatus[b] << "\t\t\t\t" << status2with3_previousEventTrigger[b] << "\t\t\t" << status2with3_previousEventNhits[b] << "\t\t\t | \t"
                << status2with3_currentEventStatus[b] << "\t\t\t" << status2with3_currentEventTrigger[b] << "\t\t" << status2with3_currentEventNhits[b] << "\t\t\t | \t"
                << status2with3_nextEventStatus[b] << "\t\t\t\t" << status2with3_nextEventTrigger[b] << "\t\t\t\t" << status2with3_nextEventNhits[b] << newline;
        }
        textOutput << newline;
        textOutput << "Details of the remaining " << counter_status2without3 << " \"Missing Trigger\" events are given below (* denotes the current event index):" << newline;
        NHitEventsMissingTrigger.push_back(counter_status2without3);
        textOutput << newline;
        textOutput << "\t" << "Status[*]" << "\t" << "Trigger[*]" << "\t" << "Nhits[*]" << newline;
        for (int b = 0; b < counter_status2without3; b++)
        {
            textOutput << "\t" << status2without3_currentEventStatus[b] << "\t\t" << status2without3_currentEventTrigger[b] << "\t\t" << status2without3_currentEventNhits[b] << newline;
        }
        textOutput << newline;

        textOutput << newline;
        textOutput << "Out of the " << counter_eventStatus3 << " \"Extra Trigger\" events, " << counter_status3without2 << " events cannot be paired with a \"Missing Trigger\" event" << newline;
        TriggerEventsMissingNhits.push_back(counter_status3without2);
        textOutput << "Details of these events are given below:" << newline;
        textOutput << newline;
        textOutput << "\t" << "Nhits" << newline;
        for (int b = 0; b < counter_status3without2; b++)
        {
            textOutput << "\t" << status3without2_currentEventNhits[b] << newline;
        }
        textOutput << newline;

        allSubruns_numberOfEvents.push_back(counter_eventStatus1);
    }

    textOutput.close();
    info << "Saved text file" << newline;

    allSubruns_nhitsVsTriggerPlots_allEvents.clear();
    allSubruns_eventStatus_allEvents.clear();
    allSubruns_eventNhits_allEvents.clear();
    allSubruns_eventTrigTypes_allEvents.clear();

    return allSubruns_numberOfEvents;
}


vector<double> DQSmellieProc::Check_3_LaserFrequency()
{
    // Prepare plots for each subrun
    vector<TH1D> allSubruns_deltaTimePlots_trigCutPassed;
    vector<TH2D> allSubruns_nhitsVsDeltaTimePlots_trigCutPassed;
    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        stringstream plot1Namestream;
        plot1Namestream << "SMELLIEDQ_Check3_deltaTimePlot_subrun" << s;
        TH1D singleSubrun_deltaTimePlot ((plot1Namestream.str()).c_str(), "Time between Successive Triggered Events for Events Passing the Trigger Cut; #delta(Time), seconds", 500, 0.0, 2.0/fSubrunLaserFrequenciesVector[s]);
        allSubruns_deltaTimePlots_trigCutPassed.push_back(singleSubrun_deltaTimePlot);

        stringstream plot2Namestream;
        plot2Namestream << "SMELLIEDQ_Check3_nhitsVsDeltaTimePlot_subrun" << s;
        TH2D singleSubrun_nhitsVsDeltaTimePlot ((plot2Namestream.str()).c_str(), "Nhits Vs. Time between Successive Triggered Events for Events Passing the Trigger Cut; Nhits; #delta(Time), seconds", 500, 0.0, 500.0, 500, 0.0, 2.0/fSubrunLaserFrequenciesVector[s]);
        allSubruns_nhitsVsDeltaTimePlots_trigCutPassed.push_back(singleSubrun_nhitsVsDeltaTimePlot);
    }
    // Calculate and plot the difference in time between successive events that passed the trigger type cut (calculate time using the 50MHz clock ticks) and fill the appropriate histograms based on the event's subrun ID
   ULong64_t previousTickCount50MHz = 0;
    if(fAllEvents_TickCounts.size() > 0){
        previousTickCount50MHz = fAllEvents_TickCounts[0];
    }
    for (unsigned int e = 1; e < fAllEvents_TrigTypes.size(); e++)
    {
        int eventSubrunID = fAllEvents_SubRunIDs[e];
        if(eventSubrunID >= fNumberOfSubruns){
            info << "Subrun number: "<<eventSubrunID<<" Greater than or eqaul  Number of subruns: "<<fNumberOfSubruns<<newline;
            continue;
        }
        int eventTriggerType = fAllEvents_TrigTypes[e];
        if (!((eventTriggerType >> fRunTriggerCutBit) & 1)) continue;

        ULong64_t currentTickCount50MHz = fAllEvents_TickCounts[e];
        ULong64_t deltaTickCount50MHz = currentTickCount50MHz - previousTickCount50MHz;
        double deltaTimeInSeconds = deltaTickCount50MHz * 2.0e-08;
        (allSubruns_deltaTimePlots_trigCutPassed[eventSubrunID]).Fill(deltaTimeInSeconds);
        (allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[eventSubrunID]).Fill(fAllEvents_Nhits[e], deltaTimeInSeconds);

        previousTickCount50MHz = currentTickCount50MHz;
    }
    // For each subrun, and if applicable, make the plots presentable and then save to the output file

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        (allSubruns_deltaTimePlots_trigCutPassed[s]).SetStats(0);
        (allSubruns_deltaTimePlots_trigCutPassed[s]).SetLineColor(kGreen+2);
        (allSubruns_deltaTimePlots_trigCutPassed[s]).SetLineWidth(2);

        (allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[s]).SetStats(0);
        (allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[s]).SetMarkerStyle(20);
        (allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[s]).SetMarkerSize(0.5);
        (allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[s]).SetMarkerColor(kGreen+2);
        WriteToRoot(&allSubruns_deltaTimePlots_trigCutPassed[s]);
        WriteToRoot(&allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[s]);

        string filename = string(allSubruns_deltaTimePlots_trigCutPassed[s].GetName());
        TCanvas outCan (filename.c_str(),filename.c_str(),1500,900);
        outCan.cd();
        allSubruns_deltaTimePlots_trigCutPassed[s].Draw();
        filename = filename+".png";
        outCan.Print(filename.c_str());

        filename = string(allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[s].GetName());
        TCanvas outCan2 (filename.c_str(),filename.c_str(),1500,900);
        outCan2.cd();
        allSubruns_nhitsVsDeltaTimePlots_trigCutPassed[s].Draw();
        filename = filename+".png";
        outCan2.Print(filename.c_str());

    }


    // For each subrun, fit a gaussian to the distribution of time differences, and find the frequency from the mean time difference
    vector<double> allSubruns_frequencies;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        TF1* deltaTimeFit = new TF1("gaus", "gaus", 0.0, 0.05);
        double deltaTime = allSubruns_deltaTimePlots_trigCutPassed[s].GetBinCenter(allSubruns_deltaTimePlots_trigCutPassed[s].GetMaximumBin());
        double laserFrequency = 1.0 / deltaTime;
        allSubruns_frequencies.push_back(laserFrequency);

        delete deltaTimeFit;
    }

    return allSubruns_frequencies;
}


vector<string> DQSmellieProc::Check_4_FibreID()
{
    // Prapare plots for each subrun
    vector<TH2D> allSubruns_firstPeak_trigCutPassed, allSubruns_secondPeak_trigCutPassed;
    vector<TH2D> allSubruns_firstPeak_trigCutPassed_pi_minus_theta, allSubruns_secondPeak_trigCutPassed_pi_minus_theta;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {

        TH2D thetaPhiPlot_firstPeak ("SMELLIEDQ_thetaPhiPlot_firstPeak", "Theta and Phi Coordinates of each Triggered PMT in the first Time Peak for Events Passing the Trigger Cut; Phi (rads); Theta (rads)", 300, -3.5, 3.5, 300, 0.0, 3.5);
        allSubruns_firstPeak_trigCutPassed.push_back(thetaPhiPlot_firstPeak);

        TH2D thetaPhiPlot_firstPeak_pi_minus_theta ("SMELLIEDQ_thetaPhiPlot_firstPeak_pi_minus_theta", "Theta and Phi Coordinates of each Triggered PMT in the first Time Peak for Events Passing the Trigger Cut; Phi (rads); #pi-Theta (rads)", 300, -3.5, 3.5, 300, 0.0, 3.5);
        allSubruns_firstPeak_trigCutPassed_pi_minus_theta.push_back(thetaPhiPlot_firstPeak_pi_minus_theta);

        TH2D thetaPhiPlot_secondPeak ("SMELLIEDQ_thetaPhiPlot_secondPeak", "Theta and Phi Coordinates of each Triggered PMT in the second Time Peak for Events Passing the Trigger Cut; Phi (rads); Theta (rads)", 300, -3.5, 3.5, 300, 0.0, 3.5);
        allSubruns_secondPeak_trigCutPassed.push_back(thetaPhiPlot_secondPeak);

        TH2D thetaPhiPlot_secondPeak_pi_minus_theta ("SMELLIEDQ_thetaPhiPlot_secondPeak_pi_minus_theta", "Theta and Phi Coordinates of each Triggered PMT in the second Time Peak for Events Passing the Trigger Cut; Phi (rads); #pi-Theta (rads)", 300, -3.5, 3.5, 300, 0.0, 3.5);
        allSubruns_secondPeak_trigCutPassed_pi_minus_theta.push_back(thetaPhiPlot_secondPeak_pi_minus_theta);
    }

    // Fill the plots of PMT locations for only events passing the trigger cut
    for (unsigned int e = 1; e < fAllEvents_TrigTypes.size(); e++)
    {
        int eventTriggerType = fAllEvents_TrigTypes[e];
        if (!((eventTriggerType >> fRunTriggerCutBit) & 1)) continue;
        int eventSubrunID = fAllEvents_SubRunIDs[e];
        if(eventSubrunID >= fNumberOfSubruns){
            info << "Subrun number: "<<eventSubrunID<<" Greater than or eqaul  Number of subruns: "<<fNumberOfSubruns<<newline;
            continue;
        }
        for (unsigned int p = 0; p < fAllEvents_pmtTimes[e].size(); p++)
        {
            if (((fAllEvents_pmtTimes[e])[p] > fFirstTimePeak_Low) && ((fAllEvents_pmtTimes[e])[p] < fFirstTimePeak_High))
            {
                (allSubruns_firstPeak_trigCutPassed[eventSubrunID]).Fill((fAllEvents_pmtPhis[e])[p], (fAllEvents_pmtThetas[e])[p]);
                (allSubruns_firstPeak_trigCutPassed_pi_minus_theta[eventSubrunID]).Fill((fAllEvents_pmtPhis[e])[p], pi-(fAllEvents_pmtThetas[e])[p]);
            }
            else if (((fAllEvents_pmtTimes[e])[p] > fSecondTimePeak_Low) && ((fAllEvents_pmtTimes[e])[p] < fSecondTimePeak_High))
            {
                (allSubruns_secondPeak_trigCutPassed[eventSubrunID]).Fill((fAllEvents_pmtPhis[e])[p], (fAllEvents_pmtThetas[e])[p]);
                (allSubruns_secondPeak_trigCutPassed_pi_minus_theta[eventSubrunID]).Fill((fAllEvents_pmtPhis[e])[p], pi-(fAllEvents_pmtThetas[e])[p]);
            }
            else continue;
        }
    }

    // For each subrun, project both the first and second peak 2D theta/phi distributions into X and Y
    vector<TH1D> allSubruns_firstPeak_thetaPlots, allSubruns_firstPeak_phiPlots, allSubruns_secondPeak_thetaPlots, allSubruns_secondPeak_phiPlots, allSubruns_firstPeak_pi_minus_thetaPlots, allSubruns_secondPeak_pi_minus_thetaPlots;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        std::stringstream projName;
        projName << "firstPeak #pi-Theta subrun"<<s;
        TH1D* pi_minus_thetaPlot_firstPeak = (allSubruns_firstPeak_trigCutPassed_pi_minus_theta[s]).ProjectionY(projName.str().c_str(),0,-1,"");
        TH1D pi_minus_thetaPlot_firstPeak_copy;
        pi_minus_thetaPlot_firstPeak->Copy(pi_minus_thetaPlot_firstPeak_copy);
        delete pi_minus_thetaPlot_firstPeak;
        pi_minus_thetaPlot_firstPeak = NULL;
        allSubruns_firstPeak_pi_minus_thetaPlots.push_back(pi_minus_thetaPlot_firstPeak_copy);

        projName.str("");
        projName.clear();

        projName << "firstPeak Theta subrun"<<s;
        TH1D* thetaPlot_firstPeak = (allSubruns_firstPeak_trigCutPassed[s]).ProjectionY(projName.str().c_str(),0,-1,"");
        TH1D thetaPlot_firstPeak_copy;
        thetaPlot_firstPeak->Copy(thetaPlot_firstPeak_copy);
        delete thetaPlot_firstPeak;
        thetaPlot_firstPeak = NULL;
        allSubruns_firstPeak_thetaPlots.push_back(thetaPlot_firstPeak_copy);

        projName.str("");
        projName.clear();

        projName << "firstPeakPhi subrun"<<s;
        TH1D* phiPlot_firstPeak = (allSubruns_firstPeak_trigCutPassed[s]).ProjectionX(projName.str().c_str(),0,-1,"");
        TH1D phiPlot_firstPeak_copy;
        phiPlot_firstPeak->Copy(phiPlot_firstPeak_copy);
        delete phiPlot_firstPeak;
        phiPlot_firstPeak = NULL;
        allSubruns_firstPeak_phiPlots.push_back(phiPlot_firstPeak_copy);

        projName.str("");
        projName.clear();
        projName << "secondPeak #pi-Theta subrun"<<s;

        TH1D* pi_minus_thetaPlot_secondPeak = (allSubruns_secondPeak_trigCutPassed_pi_minus_theta[s]).ProjectionY(projName.str().c_str(),0,-1,"");
        TH1D pi_minus_thetaPlot_secondPeak_copy;
        pi_minus_thetaPlot_secondPeak->Copy(pi_minus_thetaPlot_secondPeak_copy);
        delete pi_minus_thetaPlot_secondPeak;
        pi_minus_thetaPlot_secondPeak = NULL;
        allSubruns_secondPeak_pi_minus_thetaPlots.push_back(pi_minus_thetaPlot_secondPeak_copy);

        projName.str("");
        projName.clear();
        projName << "secondPeak Theta subrun"<<s;
        TH1D* thetaPlot_secondPeak = (allSubruns_secondPeak_trigCutPassed[s]).ProjectionY(projName.str().c_str(),0,-1,"");
        TH1D thetaPlot_secondPeak_copy;
        thetaPlot_secondPeak->Copy(thetaPlot_secondPeak_copy);
        delete thetaPlot_secondPeak;
        thetaPlot_secondPeak = NULL;
        allSubruns_secondPeak_thetaPlots.push_back(thetaPlot_secondPeak_copy);

        projName.str("");
        projName.clear();

        projName << "secondPeakPhi subrun"<<s;

        TH1D* phiPlot_secondPeak = (allSubruns_secondPeak_trigCutPassed[s]).ProjectionX(projName.str().c_str(),0,-1,"");
        TH1D phiPlot_secondPeak_copy;
        phiPlot_secondPeak->Copy(phiPlot_secondPeak_copy);
        delete phiPlot_secondPeak;
        phiPlot_secondPeak = NULL;
        allSubruns_secondPeak_phiPlots.push_back(phiPlot_secondPeak_copy);
    }

    // For each subrun, and if applicable, make the plots presentable and then save to the output file

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        (allSubruns_firstPeak_trigCutPassed[s]).SetStats(0);
        (allSubruns_firstPeak_trigCutPassed[s]).SetMarkerStyle(20);
        (allSubruns_firstPeak_trigCutPassed[s]).SetMarkerSize(0.1);
        (allSubruns_firstPeak_trigCutPassed[s]).SetMarkerColor(kGreen+2);

        (allSubruns_firstPeak_trigCutPassed_pi_minus_theta[s]).SetStats(0);
        (allSubruns_firstPeak_trigCutPassed_pi_minus_theta[s]).SetMarkerStyle(20);
        (allSubruns_firstPeak_trigCutPassed_pi_minus_theta[s]).SetMarkerSize(0.1);
        (allSubruns_firstPeak_trigCutPassed_pi_minus_theta[s]).SetMarkerColor(kGreen+2);

        (allSubruns_firstPeak_thetaPlots[s]).SetStats(0);
        (allSubruns_firstPeak_thetaPlots[s]).SetLineWidth(2);
        (allSubruns_firstPeak_thetaPlots[s]).SetLineColor(kGreen+2);
        (allSubruns_firstPeak_thetaPlots[s]).SetTitle("Theta Coordinate of each Triggered PMT in the first Time Peak for Events Passing the Trigger Cut");

        (allSubruns_firstPeak_pi_minus_thetaPlots[s]).SetStats(0);
        (allSubruns_firstPeak_pi_minus_thetaPlots[s]).SetLineWidth(2);
        (allSubruns_firstPeak_pi_minus_thetaPlots[s]).SetLineColor(kGreen+2);
        (allSubruns_firstPeak_pi_minus_thetaPlots[s]).SetTitle("#pi-Theta Coordinate of each Triggered PMT in the first Time Peak for Events Passing the Trigger Cut");

        (allSubruns_firstPeak_phiPlots[s]).SetStats(0);
        (allSubruns_firstPeak_phiPlots[s]).SetLineWidth(2);
        (allSubruns_firstPeak_phiPlots[s]).SetLineColor(kGreen+2);
        (allSubruns_firstPeak_phiPlots[s]).SetTitle("Phi Coordinate of each Triggered PMT in the first Time Peak for Events Passing the Trigger Cut");

        stringstream canvas1Namestream;
        canvas1Namestream << "SMELLIEDQ_Check4_firstPeakCanvas_subrun" << s;
        TCanvas firstPeakCanvas ((canvas1Namestream.str()).c_str(), "#pi-Theta and Phi Coordinates of each Triggered PMT in the first Time Peak for Events Passing the Trigger Cut", 1500, 900);
        firstPeakCanvas.Divide(2, 2);
        firstPeakCanvas.cd(2);
        (allSubruns_firstPeak_pi_minus_thetaPlots[s]).Draw();
        firstPeakCanvas.cd(3);
        (allSubruns_firstPeak_phiPlots[s]).Draw();

        TVirtualPad * plotPad = firstPeakCanvas.cd(4);
        plotPad->SetLogz();
        (allSubruns_firstPeak_trigCutPassed_pi_minus_theta[s]).Draw("colz");
        WriteToRoot(&firstPeakCanvas);
        string filename = string(firstPeakCanvas.GetName())+".png";
        firstPeakCanvas.Print(filename.c_str());


        (allSubruns_secondPeak_trigCutPassed[s]).SetStats(0);
        (allSubruns_secondPeak_trigCutPassed[s]).SetMarkerStyle(20);
        (allSubruns_secondPeak_trigCutPassed[s]).SetMarkerSize(0.1);
        (allSubruns_secondPeak_trigCutPassed[s]).SetMarkerColor(kGreen+2);

        (allSubruns_secondPeak_trigCutPassed_pi_minus_theta[s]).SetStats(0);
        (allSubruns_secondPeak_trigCutPassed_pi_minus_theta[s]).SetMarkerStyle(20);
        (allSubruns_secondPeak_trigCutPassed_pi_minus_theta[s]).SetMarkerSize(0.1);
        (allSubruns_secondPeak_trigCutPassed_pi_minus_theta[s]).SetMarkerColor(kGreen+2);

        (allSubruns_secondPeak_thetaPlots[s]).SetStats(0);
        (allSubruns_secondPeak_thetaPlots[s]).SetLineWidth(2);
        (allSubruns_secondPeak_thetaPlots[s]).SetLineColor(kGreen+2);
        (allSubruns_secondPeak_thetaPlots[s]).SetTitle("Theta Coordinate of each Triggered PMT in the second Time Peak for Events Passing the Trigger Cut");

        (allSubruns_secondPeak_pi_minus_thetaPlots[s]).SetStats(0);
        (allSubruns_secondPeak_pi_minus_thetaPlots[s]).SetLineWidth(2);
        (allSubruns_secondPeak_pi_minus_thetaPlots[s]).SetLineColor(kGreen+2);
        (allSubruns_secondPeak_pi_minus_thetaPlots[s]).SetTitle("#pi-Theta Coordinate of each Triggered PMT in the second Time Peak for Events Passing the Trigger Cut");

        (allSubruns_secondPeak_phiPlots[s]).SetStats(0);
        (allSubruns_secondPeak_phiPlots[s]).SetLineWidth(2);
        (allSubruns_secondPeak_phiPlots[s]).SetLineColor(kGreen+2);
        (allSubruns_secondPeak_phiPlots[s]).SetTitle("Phi Coordinate of each Triggered PMT in the second Time Peak for Events Passing the Trigger Cut");

        stringstream canvas2Namestream;
        canvas2Namestream << "SMELLIEDQ_Check4_secondPeakCanvas_subrun" << s;
        TCanvas secondPeakCanvas ((canvas2Namestream.str()).c_str(), "#pi-Theta and Phi Coordinates of each Triggered PMT in the second Time Peak for Events Passing the Trigger Cut", 1500, 900);
        secondPeakCanvas.Divide(2, 2);
        secondPeakCanvas.cd(2);
        (allSubruns_secondPeak_pi_minus_thetaPlots[s]).Draw();
        secondPeakCanvas.cd(3);
        (allSubruns_secondPeak_phiPlots[s]).Draw();

        plotPad =  secondPeakCanvas.cd(4);
        plotPad->SetLogz();
        (allSubruns_secondPeak_trigCutPassed_pi_minus_theta[s]).Draw("colz");
        WriteToRoot(&secondPeakCanvas);
        filename = string(secondPeakCanvas.GetName())+".png";
        secondPeakCanvas.Print(filename.c_str());
    }



    // For each subrun, find the Node ID using the first Timing Peak, and find the Fibre ID using the Node ID and the second Timing Peak
    vector<string> allSubruns_fibreIDs;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        double thetaOfMaxBin_firstPeak = (allSubruns_firstPeak_thetaPlots[s]).GetBinLowEdge((allSubruns_firstPeak_thetaPlots[s]).GetMaximumBin());
        double phiOfMaxBin_firstPeak = (allSubruns_firstPeak_phiPlots[s]).GetBinLowEdge((allSubruns_firstPeak_phiPlots[s]).GetMaximumBin());

        TVector3 firstPeakPosition;
        firstPeakPosition.SetMagThetaPhi(fPMTRadius, thetaOfMaxBin_firstPeak, phiOfMaxBin_firstPeak);
        double smallestDisplacement_Node = 50000.0;
        int expectedNodeID = -1;
        for (int n = 0; n < fNumberOfNodes; n++)
        {
            double currentDisplacement = (firstPeakPosition - fNodePositions[n]).Mag();
            if (currentDisplacement < smallestDisplacement_Node)
            {
                smallestDisplacement_Node = currentDisplacement;
                expectedNodeID = n;
            }
        }

        double thetaOfMaxBin_secondPeak = (allSubruns_secondPeak_thetaPlots[s]).GetBinLowEdge((allSubruns_secondPeak_thetaPlots[s]).GetMaximumBin());
        double phiOfMaxBin_secondPeak = (allSubruns_secondPeak_phiPlots[s]).GetBinLowEdge((allSubruns_secondPeak_phiPlots[s]).GetMaximumBin());

        TVector3 secondPeakPosition;
        secondPeakPosition.SetMagThetaPhi(fPMTRadius, thetaOfMaxBin_secondPeak, phiOfMaxBin_secondPeak);
        TVector3 directionFromNodeToSecondPeak = (secondPeakPosition - fNodePositions[expectedNodeID]).Unit();

        double smallestDisplacement_Fibre = 50000.0;
        int expectedFibreID = -1;
        for (int f = (3 * expectedNodeID); f < ((3 * expectedNodeID) + 3); f++)
        {
            double currentDisplacement = (directionFromNodeToSecondPeak - fFibreDirections[f]).Mag();
            if (currentDisplacement < smallestDisplacement_Fibre)
            {
                smallestDisplacement_Fibre = currentDisplacement;
                expectedFibreID = f;
            }
        }

        allSubruns_fibreIDs.push_back(fFibreIDs[expectedFibreID]);
    }

    // Clear all vectors used in this function
    allSubruns_firstPeak_trigCutPassed.clear();
    allSubruns_firstPeak_trigCutPassed_pi_minus_theta.clear();
    allSubruns_secondPeak_trigCutPassed.clear();
    allSubruns_secondPeak_trigCutPassed_pi_minus_theta.clear();
    allSubruns_firstPeak_thetaPlots.clear();
    allSubruns_firstPeak_pi_minus_thetaPlots.clear();
    allSubruns_firstPeak_phiPlots.clear();
    allSubruns_secondPeak_thetaPlots.clear();
    allSubruns_secondPeak_pi_minus_thetaPlots.clear();
    allSubruns_secondPeak_phiPlots.clear();

    return allSubruns_fibreIDs;
}


vector<double> DQSmellieProc::Check_5_LaserIntensity(vector< vector<double> > allSubruns_nhitsInformation, vector<string> allSubruns_fibreIDs)
{
    vector<double> allSubruns_intensities;

    // For each subrun ...
    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        //If a fixed wavelength Laser
        if(fLaserType[s] != "superK"){
            // ... read in the Mean Nhits vs. Intensity map for this wavelength / fibre combination from the processor's .ratdb file
            vector<double> intensitiesVector, nhitsVector;

            stringstream entryIndexStream;
            entryIndexStream << "nhits_" << allSubruns_fibreIDs[s] << "_" << fSubrunWavelengthsVector[s] ;
            string entryIndex = entryIndexStream.str();

            DBLinkPtr dbLink = DB::Get()->GetLink("DQCHECKS", "smellie");
            try
            {
                intensitiesVector = dbLink->GetDArray("intensities");
                nhitsVector = dbLink->GetDArray(entryIndex);
            }
            catch (RAT::DBNotFoundError &e) {
                info << ("DQSmellieProc::Check_5_LaserIntensity: no Intensity vs. Nhits exists for the given fibre/wavelength combination: " + entryIndex + " ... cannot calculate laser intensity Skipping subrun setting intensity to -99999");
                allSubruns_intensities.push_back(-99999);
                continue;
            }
            // ... linearly Interpolate between the values in the map to find the intensity that corresponds to this subrun's mean Nhits
            double lowIntensity = -1.0, highIntensity = -1.0, lowNhits = -1.0, highNhits = -1.0;
            for (unsigned int n = 0; n < (nhitsVector.size() - 1); n++)
            {
                if (((allSubruns_nhitsInformation[s])[1] < nhitsVector[n]) && ((allSubruns_nhitsInformation[s])[1] >= nhitsVector[n + 1]))
                {
                    lowIntensity = intensitiesVector[n];
                    highIntensity = intensitiesVector[n + 1];
                    lowNhits = nhitsVector[n];
                    highNhits = nhitsVector[n + 1];
                    break;
                }
            }

            double intensity = lowIntensity + (((highIntensity - lowIntensity) / (highNhits - lowNhits)) * ((allSubruns_nhitsInformation[s])[1] - lowNhits));
            allSubruns_intensities.push_back(intensity);
        }

        //If using the superK Laser append using the estimate of the superK laser wavelength
        else{
            // ... linearly Interpolate between the values in the map to find the intensity that corresponds to this subrun's mean Nhits
            double lowWavelength = -1.0, highWavelength = -1.0, lowNhits = -1.0, highNhits = -1.0;
            for (unsigned int n = 0; n < (fSuperKNHits.size() - 1); n++)
            {
                if (((allSubruns_nhitsInformation[s])[1] < fSuperKNHits[n]) && ((allSubruns_nhitsInformation[s])[1] >= fSuperKNHits[n + 1]))
                {
                    lowWavelength = fSuperKWavelengths[n];
                    highWavelength = fSuperKWavelengths[n + 1];
                    lowNhits = fSuperKNHits[n];
                    highNhits = fSuperKNHits[n + 1];
                    break;
                }
            }

            double superKWavelength = lowWavelength + (((highWavelength - lowWavelength) / (highNhits - lowNhits)) * ((allSubruns_nhitsInformation[s])[1] - lowNhits));
            allSubruns_intensities.push_back(superKWavelength);
        }

    }

    return allSubruns_intensities;
}


void DQSmellieProc::Check_6_DigitiserWaveforms()
{
    // Prepare a plot for each subrun
    vector<TH1D> allSubruns_singleSampleNumberPlots_trigCutPassed;

    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        stringstream plotNamestream;
        plotNamestream << "SMELLIEDQ_Check6_singleSampleNumberPlot_subrun" << s;

        TH1D singleSampleNumberPlot ((plotNamestream.str()).c_str(), "Distribution of the 100th Sample Point of each CAEN Waveform for Events Passing the Trigger Cut; Sample Value (units)", 1200, 0, 2400);
        singleSampleNumberPlot.SetLineColor(kGreen+2);
        singleSampleNumberPlot.SetLineWidth(2);
        allSubruns_singleSampleNumberPlots_trigCutPassed.push_back(singleSampleNumberPlot);
    }
    for (unsigned int e = 0; e < fAllEvents_TrigTypes.size(); e++)
    {
        int eventTriggerType = fAllEvents_TrigTypes[e];
        if (!((eventTriggerType >> fRunTriggerCutBit) & 1) ) continue;

        int eventSubrunID = fAllEvents_SubRunIDs[e];
        if(eventSubrunID >= fNumberOfSubruns){
            info << "Subrun number: "<<eventSubrunID<<" Greater than or eqaul  Number of subruns: "<<fNumberOfSubruns<<newline;
            continue;
        }
        ((allSubruns_singleSampleNumberPlots_trigCutPassed)[eventSubrunID]).Fill((fAllEvents_Waveforms[e])[fDigitiserSampleNumber]);
    }


    for (int s = 0; s < fNumberOfSubruns; s++)
    {
        WriteToRoot(&(allSubruns_singleSampleNumberPlots_trigCutPassed[s]));
        string filename = allSubruns_singleSampleNumberPlots_trigCutPassed[s].GetName();
        TCanvas outCan (filename.c_str(),filename.c_str(),1500,900);
        outCan.cd();
        allSubruns_singleSampleNumberPlots_trigCutPassed[s].Draw();
        filename = filename+".png";
        outCan.Print(filename.c_str());

    }

}