//Contact Person: Mark Stringer ms711@sussex.ac.uk

#include <TH1D.h>
#include <TVector3.h>
#include <TF1.h>
#include <TH2D.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TMath.h>

#include <RAT/DQLaserballProc.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DS/SOC.hh>
#include <RAT/DS/SOCPMT.hh>
#include <RAT/Log.hh>
#include <RAT/DS/Calib.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DS/FitResult.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/PhysicsUtil.hh>

#include <string>
#include <vector>
#include <math.h>
#include <algorithm>

using std::vector;
using std::string;

namespace RAT
{

    DQLaserballProc::DQLaserballProc() : DataQualityProc("dqlaserballProc")
    {
        //Initializing all the variables
        fDQChecks = NULL;


        fMinTime=0;
        fMaxTime=0;

        fMinQhsLow=0;
        fMaxQhsLow=0;

        fMinQhsPeak=0;
        fMaxQhsPeak=0;

        fMinQhsHigh=0;
        fMaxQhsHigh=0;

        fFitPositionThreshOnCentralAxis=0;
        fFitPositionThreshOffCentralAxis=0;

        fMinNHit=0;
        fMaxNHit=0;


        fMaxNZero=0;

        fMaxFullOccupancyMean=0;
        fMaxFullOccupancySpread=0;

        fMaxPromptOccupancyMean=0;
        fMaxPromptOccupancySpread=0;




        fMaxNBadPMTs=0;

        fMaxPMTTimeError=0;


    }

    DQLaserballProc::~DQLaserballProc()
    {
    }

    void DQLaserballProc::BeginOfRun( DS::Run& run )
    {
        //DQ begin of run
        DataQualityProc::BeginOfRun( run );
        //Getting DB
        RAT::DU::Utility::Get()->BeginOfRun();
        //loading lightpath
        fLP = RAT::DU::Utility::Get()->GetLightPathCalculator();
        //Loading GroupVelocity
        fGV = RAT::DU::Utility::Get()->GetGroupVelocity();
        //Putting PMT position database into arrays

        fPMTInfo = RAT::DU::Utility::Get()->GetPMTInfo();
        //Loading DB checks  file
        fDQChecks = DB::Get()->GetLink( "DQCHECKS", "laserball" );

        //Loading limits into variables
        fMinTime = fDQChecks->GetD("min_time");
        fMaxTime = fDQChecks->GetD("max_time");

        fMinQhsLow = fDQChecks->GetD("min_qhs_low");
        fMaxQhsLow = fDQChecks->GetD("max_qhs_low");

        fMinQhsPeak = fDQChecks->GetD("min_qhs_peak");
        fMaxQhsPeak = fDQChecks->GetD("max_qhs_peak");

        fMinQhsHigh = fDQChecks->GetD("min_qhs_high");
        fMaxQhsHigh = fDQChecks->GetD("max_qhs_high");

        fMinNHit = fDQChecks->GetI("min_nhit");
        fMaxNHit = fDQChecks->GetI("max_nhit");

        fMaxNZero = fDQChecks->GetI("max_nzero");

        fMaxFullOccupancyMean = fDQChecks->GetD("max_full_occupancy_mean");

        fMaxFullOccupancySpread = fDQChecks->GetD("max_full_occupancy_spread");

        fMaxPromptOccupancyMean = fDQChecks->GetD("max_prompt_occupancy_mean");

        fMaxPromptOccupancySpread = fDQChecks->GetD("max_prompt_occupancy_spread");


        fFitPositionThreshOnCentralAxis = fDQChecks->GetD("fit_position_thresh_on_central_axis");
        fFitPositionThreshOffCentralAxis = fDQChecks->GetD("fit_position_thresh_off_central_axis");


        fMaxNBadPMTs = fDQChecks->GetI("n_bad_pmts");

        fMaxPMTTimeError = fDQChecks->GetD("max_time_error");

        fCentralAxesX = fDQChecks->GetDArray("urm_ports_x");
        fCentralAxesY = fDQChecks->GetDArray("urm_ports_y");

        fManipInfo = DB::Get()->GetLink("CALIB_COMMON_RUN_LEVEL" ,"MANIP");

        std::vector<double> laserballPositionArray = fManipInfo->GetDArray("position");

        fSourceType = fManipInfo->GetS("source_type");
        fPositionValid = fManipInfo->GetZ("position_valid");

        fLaserballPosition = TVector3(laserballPositionArray[0],laserballPositionArray[1],laserballPositionArray[2]);
        fLaserballWavelength = fManipInfo->GetD("wavelength");



    }

    //Method to perform all checks and provide user info on whether checks failed or pass as well as setting check flags
    void DQLaserballProc::SOCAnalysis( DS::SOC& soc )
    {
        info << "DQLaserballProc::SOCAnalysis: avaliable fit names:" << newline;
        vector<string> fitNames = soc.GetFitNames();
        for( vector<string>::iterator iFitName = fitNames.begin();
                iFitName != fitNames.end(); iFitName++ )
            info << *iFitName << newline;
        info << "Now Performing Checks"<<newline;
        info << "NHitCheck"<<newline;
        int nHits = GetNHit(soc);
        AddToRATDB("mean_nhit",nHits);
        if(nHits<fMinNHit){
            info << "NHitCheck Failed too few hits NHits: "<<nHits<<" fMinNHit: "<<fMinNHit<<newline;
            Update("lbNHitCheck",false);
        }
        else if(nHits>fMaxNHit){
            info << "NHitCheck Failed too many hits NHits: "<<nHits<<" fMaxNHit: "<<fMaxNHit<<newline;
            Update("lbNHitCheck",false);
        }
        else{
            info << "NHitCheck Passed Num of NHits: "<<nHits << newline;
            Update("lbNHitCheck",true);
        }

        info << newline;

        info << "NZero Check" << newline;
        int nzero = GetNZero(soc);
        AddToRATDB("num_unhit_pmts",nzero);
        if(nzero>fMaxNZero){
            info << "NZero Check Failed: Too many unhit PMTs Unhit PMTs: "<<nzero<<" nzero_max: "<<fMaxNZero<<newline;
            Update("lbNZeroCheck",false);
        }
        else{
            info << "NZero Check Passed NZero:"<<nzero<<newline;
            Update("lbNZeroCheck",true);
        }
        info << newline;

        info << "QHS Check" << newline;
        std::vector<double>  qhsInfo = GetQHSInfo(soc);
        AddToRATDB("peak_qhs_value",qhsInfo[0]);
        if(qhsInfo[0]>fMaxQhsPeak){
            info << "QHS Test Failed: Peak value too high Peak Value: "<<qhsInfo[0]<<" fMaxQhsPeak: "<<fMaxQhsPeak<< newline;
            Update("lbqhsPeakCheck",false);
        }
        else if(qhsInfo[0]<fMinQhsPeak){
            info << "QHS Test Failed: Peak value too low Peak Value: "<<qhsInfo[0]<<" fMinQhsPeak: "<<fMinQhsPeak << newline;
            Update("lbqhsPeakCheck",false);
        }
        else{
            info << "QHS Peak Test Passed QHS Peak value: "<<qhsInfo[0]<<newline;
            Update("lbqhsPeakCheck",true);
        }
        AddToRATDB("qhs_low_half_max",qhsInfo[1]);
        if(qhsInfo[1]<fMinQhsLow){
            info << "QHS Test Failed: QHS Low too low Low Value: " <<qhsInfo[1]<<" fMinQhsLow: "<<fMinQhsLow << newline;
            Update("lbqhsLowCheck",false);
        }

        else if(qhsInfo[1]>fMaxQhsLow){
            info << "QHS Test Failed: QHS Low too high Low Value: " <<qhsInfo[1]<<" fMaxQhsLow: "<<fMaxQhsLow << newline;
            Update("lbqhsLowCheck",false);
        }
        else{
            info << "QHS Low Test Passed QHS Low Value: "<<qhsInfo[1]<<newline;
            Update("lbqhsLowCheck",true);
        }
        AddToRATDB("qhs_high_half_max",qhsInfo[2]);
        if(qhsInfo[2]<fMinQhsHigh){
            info << "QHS Test Failed: QHS high too low High Value: " <<qhsInfo[2]<<" fMinQhsHigh: "<<fMinQhsHigh << newline;
            Update("lbqhsHighCheck",false);
        }

        else if(qhsInfo[2]>fMaxQhsHigh){
            info << "QHS Test Failed: QHS high too high High Value: " <<qhsInfo[2]<<" fMaxQhsHigh: "<<fMaxQhsHigh << newline;
            Update("lbqhsHighCheck",false);
        }

        else{
            info << "QHS High Test Passed QHS High Value: "<<qhsInfo[2]<<newline;
            Update("lbqhsHighCheck",true);
        }
        info << newline;

        info << "Fit Position Test"<<newline;
        bool fitTest;
        bool onCentralAxes = false;
        //Cheking if laserball is being droped down from one of the UI ports by checking if within 20mm of one of the port positions in the x y plane
        for(unsigned int portNum=0; portNum<fCentralAxesX.size(); portNum++){
            if(TMath::Abs(fCentralAxesX[portNum]-fLaserballPosition.X())<=20  && TMath::Abs(fCentralAxesY[portNum]-fLaserballPosition.Y())<=20 ) onCentralAxes=true;
        }
        if(onCentralAxes) fitTest = FitPositionTest(soc,fFitPositionThreshOnCentralAxis);
        else fitTest = FitPositionTest(soc,fFitPositionThreshOffCentralAxis);
        if(fSourceType != "LASERBALL"){
            info << "This is not a laserball run!"<<newline;
            Update("lbPosCheck",false);

        }
        else if(!fPositionValid){
            info << "Manip position is not valid"<<newline;
            Update("lbPosCheck",false);

        }
        else if(!fitTest){
            info << "Fit Position Test Failed"<<newline;
            Update("lbPosCheck",false);
        }
        else{
            info << "Fit Position Test Passed"<<newline;
            Update("lbPosCheck",true);
        }
        info << newline;

        info << "Fit Time Test"<< newline;
        std::vector<double> timeCheck = TimeChecks(soc);
        AddToRATDB("mean_time_residual",timeCheck[0]);
        if(timeCheck[0]<-10000){
            info << "Fit Time Test Failed: Fit Failed"<<newline;
            Update("lbTimeCheck",false);
        }
        else if(timeCheck[0]<fMinTime){
            info << "Fit Time Test Failed: Time Less than Minimum Value timeCheck[0] Value: "<<timeCheck[0]<<" fMinTime: "<<fMinTime<<newline;
            Update("lbTimeCheck",false);
        }

        else if(timeCheck[0]>fMaxTime){
            info << "Fit Time Test Failed: Time More than maximum Value fit Time Value: "<<timeCheck[0]<<" fMaxTime: "<<fMaxTime<<newline;
            Update("lbTimeCheck",false);
        }
        else{
            info << "Fit Time Test Passed Time: "<<timeCheck[0]<<" ns"<<newline;
            Update("lbTimeCheck",true);
        }
        AddToRATDB("num_pmts_failing_mean_time_check",timeCheck[1]);
        if(timeCheck[1]>fMaxNBadPMTs){
            info << "Too many PMTs failing peak time Check NBadPMTs: "<<timeCheck[1]<<newline;
            Update("lbPMTMeanTimeCheck",false);
        }
        else{
            info << "Peak time Check passed  NBadPMTs: "<<timeCheck[1]<<newline;
            Update("lbPMTMeanTimeCheck",true);
        }
        AddToRATDB("num_pmts_failing_time_error_check",timeCheck[2]);
        if(timeCheck[2]>fMaxNBadPMTs){
            info << "Too many PMTs failing time Error Check NBadPMTs: "<<timeCheck[2]<<newline;
            Update("lbPMTTimeErrorCheck",false);
        }
        else{
            info << "time error Check passed  NBadPMTs: "<<timeCheck[2]<<newline;
            Update("lbPMTTimeErrorCheck",true);
        }
        info << newline;

        info << "Occupancy Test" << newline;
        std::vector<double>  occ = GetFullOccupancyMeanSpread(soc);
        Update("lbFullOccupancyMeanCheck",true);
        Update("lbFullOccupancySpreadCheck",true);
        AddToRATDB("mean_occupancy",occ[0]);
        if(occ[0]>fMaxFullOccupancyMean){
            info << "Occupancy Test Failed: Mean Full Occupancy too high Mean Occupancy: "<<occ[0]<<" fMaxOccupancyMean: "<<fMaxFullOccupancyMean<<newline;
            Update("lbFullOccupancyMeanCheck",false);
        }

        AddToRATDB("occupancy_spread",occ[0]);
        if(occ[1]>fMaxFullOccupancySpread){
            info << "Occupancy Test Failed: Occupancy Full Spread too high Occupancy: "<<occ[1]<<" fMaxOccupancySpread: "<<fMaxFullOccupancySpread<<newline;
            Update("lbFullOccupancySpreadCheck",false);
        }

        if(occ[0]<fMaxFullOccupancyMean && occ[1]<fMaxFullOccupancySpread){
            info << "Full Occupancy Test Passed"<<newline;
            info << "Full Occupancy Mean: "<<occ[0]<<newline;
            info << "FullOccupancy Spread: "<<occ[1]<<newline;
        }


        occ = GetPromptOccupancyMeanSpread(soc);
        Update("lbPromptOccupancyMeanCheck",true);
        Update("lbPromptOccupancySpreadCheck",true);
        AddToRATDB("mean_prompt_occupancy",occ[0]);
        AddToRATDB("prompt_occupancy_spread",occ[1]);
        if(occ[0]>fMaxPromptOccupancyMean){
            info << "Occupancy Test Failed: Mean Prompt Occupancy too high Mean Occupancy: "<<occ[0]<<" fMaxOccupancyMean: "<<fMaxPromptOccupancyMean<<newline;
            Update("lbPromptOccupancyMeanCheck",false);
        }

        if(occ[1]>fMaxPromptOccupancySpread){
            info << "Occupancy Test Failed: Occupancy Prompt Spread too high Occupancy: "<<occ[1]<<" fMaxOccupancySpread: "<<fMaxPromptOccupancySpread<<newline;
            Update("lbPromptOccupancySpreadCheck",false);
        }

        if(occ[0]<fMaxPromptOccupancyMean && occ[1]<fMaxPromptOccupancySpread){
            info << "Prompt Occupancy Test Passed"<<newline;
            info << "Prompt Occupancy Mean: "<<occ[0]<<newline;
            info << "PromptOccupancy Spread: "<<occ[1]<<newline;
        }


        info<<newline;
        info << "Occupancy Ratio Checks to Previous and Central run not implemented yet"<<newline;

        info << "Tests Complete"<<newline;
    }

    //End of run for processor
    void DQLaserballProc::EndOfRun( DS::Run& run )
    {
        //Getting all sources
        vector<string> SOCSourceIDs = run.GetSOCSourceIDs();
        for( vector<string>::iterator iSOCSourceID = SOCSourceIDs.begin();
                iSOCSourceID != SOCSourceIDs.end(); iSOCSourceID++ )
        {
            DS::SOC& soc = run.GetSOCData( *iSOCSourceID );
            //Performing laserball  check
            SOCAnalysis( soc );
            //DQ end of run for flags
            DataQualityProc::EndOfRun(run);
        }
    }


    //Number of PMT hits per pulse
    double  DQLaserballProc::GetNHit( DS::SOC& soc){
        double  nhits = 0;
        std::vector<unsigned int> ids = soc.GetSOCPMTIDs();

        //Iterating over all pmts summing up number of  hits
        for(unsigned int i=0; i<ids.size(); i++){
            DS::SOCPMT pmt = soc.GetSOCPMT(ids[i]);
            if(fPMTInfo.GetType(ids[i]) == DU::PMTInfo::NORMAL){
                nhits += pmt.GetTimes().size();
            }
        }
        //Divide by number of pulses to get average NHit per pulse
        return static_cast<double>(nhits)/soc.GetNPulsesTriggered();

    }

    //Number of PMTS with zero hits after entire run
    int DQLaserballProc::GetNZero(DS::SOC& soc){
        int nzero = 0;
        //Iterating through all PMTs to get number of unhit pmts
        for(unsigned int i=0; i<fPMTInfo.GetCount(); i++){
            //If PMT is normal and doesnt exist in socfile the pmt must be unhit
            if(fPMTInfo.GetType(i) == DU::PMTInfo::NORMAL && !soc.SOCPMTExists(i)){
                nzero++;
            }
        }
        return nzero;
    }

    //Method to check the Mean occupancy and the spread on it
    //Full Occupancy (number of Hits on a pmt)/number of pulses
    std::vector<double> DQLaserballProc::GetFullOccupancyMeanSpread(DS::SOC& soc ){
        //Mean and spread
        std::vector<double>output(2,0);
        //Setting up occupancy histo
        TH1D  occupancy("LBDQ_Full_Occupancy_of_the_PMTs","LBDQ_Full_Occupancy_of_the_PMTS",50,0,0.3);
        //Filling histogram with occupancy
        for(unsigned int i=0; i<fPMTInfo.GetCount(); i++){
            if(fPMTInfo.GetType(i) == DU::PMTInfo::NORMAL){
                if(soc.SOCPMTExists(i)){
                    DS::SOCPMT pmt = soc.GetSOCPMT(i);
                    occupancy.Fill(((double)pmt.GetTimes().size())/soc.GetNPulsesTriggered());
                }
                //If pmt does not exist it is unhit hence fill histogram with 0
                else{
                    occupancy.Fill(0);
                }
            }
        }
        //Returning mean value and RMS of histogram
        output[0] = occupancy.GetMean();
        output[1] = occupancy.GetRMS();
        WriteToRoot(&occupancy);
        return output;
    }

    //Method to get prompt occupancy
    std::vector<double> DQLaserballProc::GetPromptOccupancyMeanSpread(DS::SOC& soc ){
        //Mean and spread
        std::vector<double>output(2,0);
        //Setting up occupancy histo
        TH1D  occupancy ("LBDQ_Prompt_Occupancy_of_the_PMTs","LBDQ_Prompt_Occupancy_of_the_PMTS",50,0,0.3);
        //Filling histogram with occupancy
        for(unsigned int i=0; i<fPMTInfo.GetCount(); i++){
            if(fPMTInfo.GetType(i) == DU::PMTInfo::NORMAL){
                if(soc.SOCPMTExists(i)){
                    DS::SOCPMT pmt = soc.GetSOCPMT(i);
                    occupancy.Fill((double)pmt.GetPromptOccupancy()/soc.GetNPulsesTriggered());
                }
                //If pmt does not exist it is unhit hence fill histogram with 0
                else{
                    occupancy.Fill(0);
                }
            }
        }
        //Returning mean value and RMS of histogram
        output[0] = occupancy.GetMean();
        output[1] = occupancy.GetRMS();
        WriteToRoot(&occupancy);
        return output;
    }



    //FitPosition check
    bool DQLaserballProc::FitPositionTest(DS::SOC& soc, double threshold){
        //Get Calib data
        //How many fits are we going to perform?
        vector<string> fitNames = soc.GetFitNames();
        TVector3 fitpos;
        for( vector<string>::iterator iFitName = fitNames.begin(); iFitName != fitNames.end(); iFitName++ ){
            DS::FitResult fit = soc.GetFitResult(*iFitName);
            info << "Fit Position threshold is: "<<threshold<<newline;
            //Checking if fit is valid
            if(!fit.GetValid()){
                info << "Fit result is invalid DQ will fail" << newline;
                return false;
            }

            //Should only be one vertex, the laserball source position (CHECK THIS?)
            if(fit.GetVertexCount()>1){
                info << "Should only be one vertex"<<newline;
                return false;
            }


            //Now getting fit vertex
            DS::FitVertex vert = fit.GetVertex(0);
            //Checking if position is valid
            if(!vert.ValidPosition()){
                info << "Fit Position is not Valid" << newline;
                return false;
            }
            //Getting position of fit
            fitpos = vert.GetPosition();
            //Calculating difference  between calib position and fit position
            TVector3 diff = fLaserballPosition-fitpos;
            //Checking if Test Fails
            if(diff.Mag()>threshold){
                info << " Fit vertex "<< *iFitName << " not close enough to laserball position : FAIL." << newline;
                info << "Fit Position is: ("<<fitpos.X()<<","<<fitpos.Y()<<","<<fitpos.Z()<<")"<<newline;
                info << "Calib Position is: ("<<fLaserballPosition.X()<<","<<fLaserballPosition.Y()<<","<<fLaserballPosition.Z()<<")"<<newline;
                return false;
            }

        }
        //Test passed
        info  << "Fit Position test passed." << newline;
        info << "Fit Position is: ("<<fitpos.X()<<","<<fitpos.Y()<<","<<fitpos.Z()<<")"<<newline;
        info << "Manip Position is: ("<<fLaserballPosition.X()<<","<<fLaserballPosition.Y()<<","<<fLaserballPosition.Z()<<")"<<newline;
        //Add to RATDB
        std::vector<double> manipPos;
        manipPos.push_back(fLaserballPosition.X());
        manipPos.push_back(fLaserballPosition.Y());
        manipPos.push_back(fLaserballPosition.Z());
        std::vector<double> fitPosition;
        fitPosition.push_back(fitpos.X());
        fitPosition.push_back(fitpos.Y());
        fitPosition.push_back(fitpos.Z());
        AddToRATDB("manip_laserball_position",manipPos);
        AddToRATDB("fit_laserball_position",fitPosition);
        return true;

    }






    //Method to measure time difference between hit time (Shifted for all electronics and TOF) and firing time
    //returns
    //>Mean Time offset over all pmts
    //>Number of PMTS which fail as their peak time falls outside limits
    //>Number of PMTS which fail as their peak error is too large
    std::vector<double> DQLaserballProc::TimeChecks(DS::SOC& soc){
        double meanTimeFails = 0;
        double errorTimeFails = 0;
        //Creating histogram to store peak prompt values
        TH1D  timeHisto  ("LBDQ_TimeHisto","LBDQ_TimeHisto",125,-20,20);
        TH1D  pmtTime ("LBDQ_PMT_Mean_Time_Check","LBDQ_PMT_Mean_Time_Check",10000,0,10000);
        TH1D  pmtTimeError  ("LBDQ_PMT_Error_Time_Check","LBDQ_PMT_Error_Time_Check",10000,0,10000);
        //Getting time offset due to electronics
        double globalTimeOffset = soc.GetGlobalTimeOffset();
        //getting wavelength of lb and converting wavelength into mm
        double wavelength = fLaserballWavelength*1.0e-6;
        //Getting energy of photons
        double energy = RAT::util::WavelengthToEnergy(wavelength);
        //Iterating through PMTs
        for(unsigned int i=0; i<fPMTInfo.GetCount(); i++){
            //Checking if PMT exists
            if(soc.SOCPMTExists(i) ){
                //Loading PMT and getting position from database
                DS::SOCPMT pmt = soc.GetSOCPMT(i);
                //Checking if a regular PMT
                if(fPMTInfo.GetType( pmt.GetID() ) == DU::PMTInfo::NORMAL && pmt.GetPeakFindOK()==DS::SOCPMT::eSuccess  ){

                    int lcn = pmt.GetLCN();
                    TVector3 pmtPos = fPMTInfo.GetPosition(pmt.GetID());
                    TVector3 pmtDir = fPMTInfo.GetDirection(pmt.GetID());
                    //Lightpath calculator calculating with tolerance of 10mm
                    fLP.CalcByPosition(fLaserballPosition,pmtPos,energy,10);
                    //Getting distances travelled in detector
                    double distInWater = fLP.GetDistInWater();
                    double distInScint = fLP.GetDistInInnerAV();
                    double distInAV = fLP.GetDistInAV();
                    double angleOfEntry = acos(fLP.GetIncidentVecOnPMT().Dot(pmtDir));
                    angleOfEntry = angleOfEntry*TMath::RadToDeg();
                    //Using Group Vel to get time of flight
                    double tof = fGV.CalcByDistance(distInScint,distInAV,distInWater,energy)+fGV.PMTBucketTime(angleOfEntry);
                    //Hit time - time of flight - electronics offset
                    char name[17];
                    sprintf(name,"LBDQ_PMT%d",lcn);
                    TH1D pmtHisto(name,name,10000,-1000,1000);
                    double meanTime = pmt.GetTimeCentroid()-tof-globalTimeOffset;
                    double timeError = pmt.GetTimeCentroidError();
                    if(meanTime>fMaxTime || meanTime<fMinTime){
                        pmtTime.SetBinContent(lcn,1.0);
                        meanTimeFails++;

                    }
                    if(timeError>fMaxPMTTimeError){
                        pmtTimeError.SetBinContent(lcn,1.0);
                        errorTimeFails++;
                    }
                    vector<float> tacs = pmt.GetTimes();
                    for(unsigned int k=0; k<tacs.size(); k++){
                        float TACtime = tacs[k]-tof-globalTimeOffset;
                        if(TACtime>-1000){
                            pmtHisto.Fill(TACtime);
                            timeHisto.Fill(TACtime);
                        }
                    }

                }
            }
        }
        //Drawing output plot and fitting to gaussian

        TCanvas  c2  ("LBDQ_TimeHistogram"); timeHisto.Fit("gaus");
        timeHisto.Draw();
        WriteToRoot(&c2);
        TCanvas  c3  ("LBDQ_pmtTimeCrateMap");
        TH2D pmtTimeCrateMap = CrateMapFromChannelHisto(pmtTime);
        pmtTimeCrateMap.Draw("colz");
        WriteToRoot(&c3);
        TH2D pmtTimeErrorCrateMap = CrateMapFromChannelHisto(pmtTimeError);
        TCanvas c4 ("LBDQ_pmtTimeErrorCrateMap");
        pmtTimeErrorCrateMap.Draw("colz");
        WriteToRoot(&c4);
        TF1 * func = timeHisto.GetFunction("gaus");
        //Returing mean value of fitted gaussian
        std::vector<double> output (3,0);
        output[0] = -10000;
        if(func!=NULL){
            output[0] = func->GetParameter(1);
        }
        output[1] = meanTimeFails;
        output[2] = errorTimeFails;
        delete func;
        return output;

    }


    //Method to obtain
    //>The peak value of QHS histogram (QHS Peak)
    //>The half maximum value on the left side of the peak (QHS Low)
    //>The half maximum value on the right side of the peak (QHS High)
    std::vector<double> DQLaserballProc::GetQHSInfo(DS::SOC& soc){
        //Loading the QHS Histogram
        TH1D qhsHisto ("LBDQ_QHS_Histogram","LBDQ_QHS_Histogram",500,0,500);
        //Filling QHS Histogram
        for(unsigned int i=0; i<fPMTInfo.GetCount(); i++ ){
            if(soc.SOCPMTExists(i)){
                DS::SOCPMT pmt = soc.GetSOCPMT(i);
                std::vector<float> qhsVals = pmt.GetQHSs();
                for(unsigned int j=0; j<qhsVals.size(); j++){
                    qhsHisto.Fill(qhsVals[j]);
                }
            }
        }
        //Drawing the QHS histogram to output
        TCanvas  c1 ("LBDQ_QHSHistogram");
        qhsHisto.Draw();
        //Now performing analysis on histogram
        //Finding peak value by fitting to gaussian
        //Using limits which seem to work well with simulation data may need to change
        qhsHisto.Fit("gaus","","",15,35);
        //Drawing function on plot
        TF1 * peakFunc = qhsHisto.GetFunction("gaus");
        peakFunc->SetLineColor(1);
        peakFunc->Draw("same");
        //Getting height of peak (norm) and peak value (mean)
        double mean = peakFunc->GetParameter(1);
        double norm = peakFunc->GetParameter(0);
        int peakBin = qhsHisto.GetXaxis()->FindFixBin(mean);
        int lowBin = -1;
        int upBin = 999999;
        //Getting 25% and 75% bins on left side of peak
        for(int i=1; i<peakBin; i++){
            if(qhsHisto.GetBinContent(i)>(0.75*norm) && upBin > peakBin) upBin = i-1;
            if(lowBin < 0 && qhsHisto.GetBinContent(i)>(0.25*norm)) lowBin = i;
        }
        //Setting up linear function between the 25% and 75% bins to interpolate and get the 50% value
        TF1  lowerFunc ("lowerFunc","([1]*x)+[0]",qhsHisto.GetBinCenter(lowBin),qhsHisto.GetBinCenter(upBin));
        lowerFunc.SetLineColor(2);
        //Fitting function
        qhsHisto.Fit("lowerFunc","+","",qhsHisto.GetBinCenter(lowBin),qhsHisto.GetBinCenter(upBin));
        lowerFunc.Draw("same");
        //Getting gradient and intercept of function
        double lowIntercept = lowerFunc.GetParameter(0);
        double lowGrad = lowerFunc.GetParameter(1);
        //Rearranging equation to solve for 50% value
        double qhsLow = (0.5*norm-lowIntercept)/lowGrad;

        //Now doing the same for the right side of the peak
        lowBin=-1;
        upBin = -1;
        for(int i=peakBin; i<qhsHisto.GetXaxis()->GetNbins(); i++){
            if(qhsHisto.GetBinContent(i)<(0.75*norm) && lowBin==-1) lowBin = i-1;
            if(upBin==-1 && qhsHisto.GetBinContent(i)<(0.25*norm)) upBin = i-1;
        }
        TF1  upperFunc ("upperFunc","([1]*x)+[0]",qhsHisto.GetBinCenter(lowBin),qhsHisto.GetBinCenter(upBin));
        qhsHisto.Fit("upperFunc","+","",qhsHisto.GetBinCenter(lowBin),qhsHisto.GetBinCenter(upBin));
        upperFunc.SetLineColor(2);
        upperFunc.Draw("same");
        double upperIntercept = upperFunc.GetParameter(0);
        double upperGrad = upperFunc.GetParameter(1);
        double qhsHigh = (0.5*norm-upperIntercept)/upperGrad;
        std::vector<double>output(3,0);
        //Outputting peak then low then high
        output[0] = mean;
        output[1] = qhsLow;
        output[2] = qhsHigh;
        WriteToRoot(&c1);
        delete peakFunc;
        return output;
    }

    //If channel fails input histobin is set to 1 -> Red In output crate map
    //If channel passed input histobin -> Blue in output crate map
    TH2D DQLaserballProc::CrateMapFromChannelHisto(TH1D channelFails){
        char histoName [100];
        sprintf(histoName,"LBDQ_%s_Crate_Map",channelFails.GetName());
        TH2D outputHisto(histoName,histoName,280,0,1,100,0,1);
        int crateGap = 10;
        int xOffset = 10;
        int yOffset = 10;

        for(unsigned int i=0; i<fPMTInfo.GetCount(); i++){
            int crateNum = crateNumber(i);
            int cardNum = cardNumber(i);
            int channelNum = channelNumber(i);

            int ybin = channelNum+yOffset;
            int xbin = cardNum+(crateNum*(16+crateGap))+xOffset;
            if( crateNum >= 9){
                ybin += 50;
                xbin = cardNum+((crateNum-9)*(16+crateGap))+xOffset;
            }
            //16 cards per crate
            int bin = outputHisto.GetBin(xbin,ybin);
            outputHisto.SetBinContent(bin,channelFails.GetBinContent(i)*10000+1);
        }
        outputHisto.SetStats(false);
        return outputHisto;
    }

    int DQLaserballProc::crateNumber(int lcn){
        return lcn/512;
    }

    int DQLaserballProc::cardNumber(int lcn){
        return ((lcn-(512*crateNumber(lcn)))/32);
    }

    int  DQLaserballProc::channelNumber(int lcn){
        return lcn-(32*cardNumber(lcn))-(512*crateNumber(lcn));
    }
} //namespace RAT