#include "TCanvas.h"
#include "TChain.h"
#include "TGraphErrors.h"
#include "TAxis.h"
#include "TCutG.h"
#include "TH2F.h"
#include "TMath.h"
#include <utility>
#include <math.h>

using namespace std;

void SignalToNoiseRatio();

void GetCountPerThreshold(Int_t* eventCountPerThreshold, Float_t* cuts, const Int_t numCuts,  const Float_t DCVEne);

Float_t eventEnergy(const Float_t *detectorChannelArrayEnergy, const Int_t numChannelsInModule, const Int_t startIndex = 0);

int main(){
    SignalToNoiseRatio();

    return 0;
}

void SignalToNoiseRatio(){
    const Int_t NUM_CUTS = 21;
    Float_t DCVModuleEne[32] = {};
    Float_t DCVEne = 0;
    Float_t Pi0Pt = 0;
    Float_t Pi0RecZ = 0;
    Float_t currentRatio = 0;

    Int_t countPerThresholdKLpi0nunu[NUM_CUTS] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    Int_t countPerThresholdKLpipipi0[NUM_CUTS] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    
    Float_t xCoordSignal[] = {3000, 4000, 4700, 4700, 3000, 3000};
    Float_t yCoordSignal[] = {130, 130, 150, 250, 250, 130};

    Float_t xVals[NUM_CUTS] = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
                               1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0};
    Float_t yVals[NUM_CUTS] = {0};

    Float_t xErrors[NUM_CUTS];
    Float_t yErrors[NUM_CUTS];

    TCanvas* c = new TCanvas("c", "c", 600, 600);

    TChain* chainKLpipipi0 = new TChain("RecTree");
    TChain* chainKLpi0nunu = new TChain("RecTree");

    TGraph* StoN = new TGraph(NUM_CUTS, xVals, yVals);

    chainKLpipipi0->Add("../MCfiles/RUN82_KLpipipi0_g2ana_KinCuts_VetoCuts_NoDCV_NoBPCV_NoBHCV_NoCC06_NoCC05_NoCC04.root");
    chainKLpi0nunu->Add("../MCfiles/RUN82_KLpi0nunu_g2ana_KinCutsNoPi0PtNoPi0RecZ_NoDownStreamVeto.root");

    chainKLpipipi0->SetBranchAddress("DCVModuleEne", DCVModuleEne);
    chainKLpipipi0->SetBranchAddress("Pi0Pt", &Pi0Pt);
    chainKLpipipi0->SetBranchAddress("Pi0RecZ", &Pi0RecZ);

    printf("Reading KLpipipi0 chain\n");

    for(Int_t i = 0; i < chainKLpipipi0->GetEntries(); i++){
        chainKLpipipi0->GetEntry(i);

        DCVEne = eventEnergy(DCVModuleEne, 32);

        GetCountPerThreshold(countPerThresholdKLpipipi0, xVals, NUM_CUTS, DCVEne);

    }

    chainKLpi0nunu->SetBranchAddress("DCVModuleEne", DCVModuleEne);
    chainKLpi0nunu->SetBranchAddress("Pi0Pt", &Pi0Pt);
    chainKLpi0nunu->SetBranchAddress("Pi0RecZ", &Pi0RecZ);

    printf("Reading KLpi0nunu chain\n");

    for(Int_t i = 0; i < chainKLpi0nunu->GetEntries(); i++){
        chainKLpi0nunu->GetEntry(i);

        DCVEne = eventEnergy(DCVModuleEne, 32);

        GetCountPerThreshold(countPerThresholdKLpi0nunu, xVals, NUM_CUTS, DCVEne);

    }

    // Construct the S/N value
    for(Int_t i = 0; i < StoN->GetN(); i++){
        if (countPerThresholdKLpi0nunu[i] == 0 || countPerThresholdKLpi0nunu[i] + countPerThresholdKLpipipi0[i] == 0){
            currentRatio = 0;
        }
        else {
            currentRatio = countPerThresholdKLpi0nunu[i] / sqrt(countPerThresholdKLpi0nunu[i] + countPerThresholdKLpipipi0[i]);   
        }

        printf("i: %i\n", i);
        printf("cut val: %f\n", xVals[i]);
        printf("ratio val: %f\n", currentRatio);

        StoN->SetPoint(i, (Double_t)xVals[i], (Double_t)currentRatio);
    }

    StoN->SetFillColorAlpha(6, 0.5);

    StoN->SetTitle("Signal to Noise Ratio Of Each Threshold");

    StoN->GetXaxis()->SetTitle("Energy Threshold [MeV]");

    StoN->SetMarkerStyle(8);

    StoN->SetMarkerColorAlpha(6, 0.5);

    StoN->Draw("PA");

    c->Update();

    c->Print("../DCVcuts/plots/StoN-plot.pdf");

    delete c;
    delete chainKLpi0nunu;
    delete chainKLpipipi0;
}

void GetCountPerThreshold(Int_t eventCountPerThreshold[], Float_t cuts[], const Int_t numCuts,  const Float_t DCVEne){
    for(Int_t i = 0; i < numCuts; i++){
        if(DCVEne < cuts[i]){
            eventCountPerThreshold[i]++;
        }
    }
}

Float_t eventEnergy(const Float_t *detectorChannelArrayEnergy, const Int_t numChannelsInModule, const Int_t startIndex){
    Float_t sum = 0;

    for(Int_t i = startIndex; i < numChannelsInModule + startIndex; i++){
        sum += detectorChannelArrayEnergy[i];
    }

    return sum;
}