////////////////////////////////////////////////////////////////////
// Plot the event energies for the different types of events.
#include "readEvents.C"

void plotEnergy() {
    readEvents();
    TH1F* dataEnergy = new TH1F("dataEnergy","Energy of Two Ring Events",
                                40,0.0,1335.0);
    dataEnergy->SetXTitle("Visible Energy (MeV)");

    TH1F* mcEnergy = new TH1F("mcEnergy","Energy of Two Ring Events",
                                40,0.0,1335.0);
    mcEnergy->SetXTitle("Visible Energy (MeV)");

    TH1F* pi0Energy = new TH1F("pi0Energy","Energy of Two Ring Events",
                                40,0.0,1335.0);
    pi0Energy->SetXTitle("Visible Energy (MeV)");

    for (struct eventData* ev = gEventData; 
         ev < &gEventData[gNEvents]; 
         ++ev) {
        if (ev->type == 0) {
            dataEnergy->Fill(ev->energy);
        }
        else if (ev->type == 1) {
            mcEnergy->Fill(ev->energy);
        }
        else if (ev->type == 2) {
            mcEnergy->Fill(ev->energy);
            pi0Energy->Fill(ev->energy);
        }
    }

    dataEnergy->SetMinimum(0);
    dataEnergy->Draw();
    gPad->Print("energy-data.png");

    mcEnergy->SetMinimum(0);
    mcEnergy->Draw();
    pi0Energy->SetLineStyle(2);
    pi0Energy->SetLineColor(2);
    pi0Energy->Draw("same");
    TLegend *mcLegend = new TLegend(0.4,0.7,0.65,0.85);
    mcLegend->AddEntry(mcEnergy,"All MC Events","L");
    mcLegend->AddEntry(pi0Energy,"#pi^{0} MC Events","L");
    mcLegend->Draw();

    gPad->Print("energy-mc.png");

    mcEnergy->Scale(1490.0/36500);
    mcEnergy->SetLineStyle(2);
    dataEnergy->Draw("E");
    mcEnergy->Draw("same");
    TLegend *overlay = new TLegend(0.4,0.7,0.65,0.85);
    overlay->AddEntry(dataEnergy,"Data Events","L");
    overlay->AddEntry(mcEnergy,"MC Events","L");
    overlay->Draw();
    gPad->Print("energy-overlay.png");

    double chi2 = 0;
    for (int i=1; i<=dataEnergy->GetNbinsX(); ++i) {
        double data = dataEnergy->GetBinContent(i);
        double mc = mcEnergy->GetBinContent(i);
        double diff = data-mc;
        chi2 += diff*diff/mc;
    }
    std::cout << "Pearson's " << chi2 << std::endl;
    std::cout << "Pearson's Probability " << TMath::Prob(chi2,dataEnergy->GetNbinsX()) << std::endl;
    

}