#include <vector>
#include <cstring>
#include <iostream>
#include <fstream>
//#include "/usr/include/root/TH2.h"
#include "TH2.h"
#include <TH2.h>
#include "TCanvas.h"
#include "TFile.h"
#include <math.h>
#include "TLegend.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TLine.h"
#include "TLatex.h"
#include "TStyle.h"
//#define draw_strips main

//this requires ROOT, I generally just run root draw_strips.cc+, but there are probably better options
//This prints histograms of the covariance map based on number of hits, the expected covariance map based on diagonal values, and a map of how many stdevs observed is away from expected
//This also prints 1d histogram distributions of said number of stdevs between expected and observed  
void draw_strips()
//int main(void)
{
    // typical making a canvas and setting the color scheme, palette can be switched. A picture of what this palette makes is attached in the email
    TCanvas *c = new TCanvas("pad1","pad1",150,10,990,660);;
    gStyle->SetPalette(55);
    //FIXME below is the file name of the csv made in covariance.cc, make sure it matches up. I just replace and recompile because I can't get inputs with argv to work with ROOT. You're welcome to try though
    std::ifstream in("/home/rutca/c++/test/covariance_dac11_new_big.csv"); 
    // the following are 4 2d histograms I make for each file, 240 by 240
    
    TH2F * output = new TH2F("covariance","covariance",240, 0.5,240.5, 240, 0.5, 240.5);
    // output is the data collected
    
    TH2F * expected = new TH2F("expected","expected",240, 0.5,240.5, 240, 0.5, 240.5);
    // expected is the extrapolation from the diagonals 
    // e.g the number of hits at 1,2 should be the number at 1,1 * 2,2 divided by (total events collected * 4). The 4 is because diagonals are double counted
    
    TH2F * difference = new TH2F("difference","difference",240, 0.5,240.5, 240, 0.5, 240.5);
    // difference is the number of standard deviations between expected and experimental (poisson distribution so stdev is sqrt of output)
    // Be careful with high dac values, since some positions end with very low chances of having a hit, so if they get 1, it ends up being hundreds of deviations out from expected
    
    TH2F * strange = new TH2F("outliers","outliers",240, 0.5,240.5, 240, 0.5, 240.5);
    // this is a map of points with high difference values, so this would in theory map positions with high or low covariance


    gStyle->SetPalette(55); //? not sure why this is here, could probably be removed

    // the following are 1d 
    TH1F * ssa1 = new TH1F("ssa1 distribution","ssa1 distribution",100,-10,10);
    // This is a histogram of the difference (the 2d hist) for the first ssa. This should look like a gaussian if things go right


    TH1F * ssa2 = new TH1F("ssa2 distribution","ssa2 distribution",100,-10,10);
    // Same as above but for ssa2


    TH1F * overlap = new TH1F("overlap distribution","overlap distribution",100,-10,10);
    // same as above but for the overlap region 

    TH1F * overall = new TH1F("overall distribution","overall distribution",100,-10,10);
    // same as above but for the entire 240 x 240 difference histogram

    // getting rid of stat boxes
    output->SetBit(TH1::kNoStats);
    ssa1->SetBit(TH1::kNoStats);
    ssa2->SetBit(TH1::kNoStats);
    overlap->SetBit(TH1::kNoStats);
    overall->SetBit(TH1::kNoStats);
//    expected->SetBit(TH1::kNoStats);
    difference->SetBit(TH2::kNoStats);

    //initializing to read the input csv
    std::string line;
    int comma;
    std::string temp;
    int hits;
    int x = 0;
    int y = 0;
    std::vector<float> diagonals;
    // below fills in output and a vector of diagonals for use in expected
    for(int j = 1; j <= 240; j++)
    {
        getline(in, line);
        for(int i = 1; i <=240; i++)
        {

            comma = line.find(",");
            temp = line.substr(0,comma);
//            std::cout << temp << ",";
            hits = std::stoi(temp);
            if(i == j)
            {
                diagonals.push_back(hits);
            }
            if(hits!=0)
            {
                output->Fill(j,i,hits);
            }
            line.erase(0,comma+1);
        }
//        std::cout << std::endl;
    }
    for(int i = 0; i <240; i++) // debug purposes 
    {
//        std::cout <<output->GetBinContent(120,i) << ",";
    }
    std::cout << std::endl;

    // below fills in expected, difference and outlier
    for(int i = 1; i <=240; i ++)
    {
        for(int j = 1; j<=240;j++)
        {

            expected->Fill(i,j,diagonals[i-1]*diagonals[j-1]/20000000); // 20,000,000 since 5,000,000 events, and each diagonal is counted twice
            if(abs(i-j) > 1)
            {
                float diff = output->GetBinContent(i,j)-expected->GetBinContent(i,j);
                float dev = sqrt(output->GetBinContent(i,j));
                difference->Fill(i,j,diff/dev);

                if(i == 40 && j == 45) // testing purposes, see if things are being calculated properly or why a position has weird behavior
                {
                    std::cout << "diff:" << diff << std::endl;
                    std::cout << "dev:" << dev << std::endl;
                    std::cout << "expected:" << expected->GetBinContent(i,j)<< std::endl;
                    std::cout << "observed:" << output->GetBinContent(i,j) << std::endl;
                    std::cout << "xvalue:" << diagonals[i-1]<<std::endl;
                    std::cout << "yvalue:" << diagonals[j-1] << std::endl;
                }
                if(abs(diff/dev) > 10) // 10 stdev away, outside the range of the 1d distribution histograms
                {
                    strange->Fill(i,j);
                }
                
                overall->Fill(difference->GetBinContent(i,j)); // filling overall
//                std::cout << difference->GetBinContent(i,j) << ",";
            }
        }
//        std::cout << std::endl;
    }

    // below fills the 1d distribution histograms
    for(int i = 1; i <=120; i++)
    {
        for(int j = 1; j <=120; j++)
        {
            ssa1->Fill(difference->GetBinContent(i,j));
            overlap->Fill(difference->GetBinContent(i,j+120));
            ssa2->Fill(difference->GetBinContent(i+120,j+120));
        }
    }
    // drawing all the relevant pieces 
    difference->Draw("COLZ");
    TCanvas *c2 = new TCanvas("pad2","pad2",150,10,990,660);
    ssa1->Draw("HIST");
    TCanvas *c3 = new TCanvas("pad3","pad3",150,10,990,660);
    ssa2->Draw("HIST");
    TCanvas *c4 = new TCanvas("pad4","pad4",150,10,990,660);
    overlap->Draw("HIST");
    TCanvas *c5 = new TCanvas("pad5","pad5",150,10,990,660);
    overall->Draw("HIST");
    TCanvas *c6 = new TCanvas("pad6","pad6",150,10,990,660);
    strange->Draw("COLZ");

}