//
//  actual.c
//
//
//  Created by Damien Storimans on 13/06/16.
//
//

#include "TInterpreter.h"
#include "TCanvas.h"
#include "TSystem.h"
#include "TFile.h"
#include "TH2.h"
#include "TNtuple.h"
#include "TPaveLabel.h"
#include "TPaveText.h"
#include "TFrame.h"
#include "TSystem.h"
#include "TInterpreter.h"
#include "TMath.h"

#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <iomanip>
#include <cmath>

using std::cout;
using std::vector;
using std::ofstream;
using std::string;
using std::setprecision;
using std::min_element;
#include <stdio.h>
//get branches for int(from claire)

vector<Int_t> GetBranchInt ( TString rootfile ,
                            TString branchname )
{
    TFile f(rootfile.Data());
    TTree *T = (TTree*)f.Get("particle");
    Long_t n_entries = (Long_t)T->GetEntries();
    vector<Int_t> branch(n_entries);
    Int_t temp;
    TBranch *branch_p = T->GetBranch(branchname);
    branch_p-> SetAddress (&temp);
    for (Long_t i=0; i<n_entries ; ++i)
    {
        branch_p->GetEntry(i);
        branch[i] = temp;
    }
    f.Close ();
    return branch ;
}
//get branch for double (from claire)
vector < Double_t > GetBranchDouble ( TString rootfile ,
                                     TString branchname )
{
    TFile f(rootfile.Data());
    TTree *T = (TTree*)f.Get("particle");
    Long_t n_entries = (Long_t)T->GetEntries();
    vector<Double_t> branch(n_entries);
    Double_t temp;
    TBranch *branch_p = T->GetBranch(branchname);
    branch_p-> SetAddress (&temp);
    for (Long_t i =0; i <n_entries ; ++i)
    {
        branch_p->GetEntry(i);
        branch[i] = temp;
    }
    f.Close ();
    return branch ;
}//find total energy (adapted from claire)
vector<Double_t> TotalE ( string file )
{
    
    vector<Int_t> type = GetBranchInt(file,"ParticleID");
    vector<Double_t> px = GetBranchDouble(file,"Px");
    vector<Double_t> py = GetBranchDouble(file,"Py");
    vector<Double_t> pz = GetBranchDouble(file,"Pz");
    

    //get stuff from tree

    
    vector<Double_t> E;
    E.clear();


    //compute mc+pc
    Double_t momentum_sq;
    momentum_sq=0.;
    Double_t energy;
    energy=0.;
    int sizetype = type.size();
    for (long i=0; i<sizetype; ++i)
    {
        //SELECT MUON:
        if (type[i]==5 || type[i]==6) // Muon + or -
        {
        // OR DONT -> Does this makes sense without selecting a particle?
        // muon (mc)^2: 0.10565837**2
            // Energy in GeV : mass GeV / c ^2 and momentum GeV / c
            momentum_sq = std::pow(px[i],2)+ std::pow(py[i],2)+ std::pow(pz[i],2);
            energy = TMath::Sqrt(std::pow(0.1056583715,2)+momentum_sq );
            E.push_back(energy);
        }
    }
    return E;
}
vector<Double_t> XYenergy( string file)
{
    //get stuff from tree
    vector<Int_t> type = GetBranchInt(file,"ParticleID");
    vector<Double_t> px = GetBranchDouble(file,"Px");
    vector<Double_t> py = GetBranchDouble(file,"Py");
    vector<Double_t> pz = GetBranchDouble(file,"Pz");
    vector<Double_t> Et;
    
    //computes sqrt(px^2+py^2)
    
    Double_t momentum_xy, energyxy;
    int sizetype = type.size();
    for (long i=0; i<sizetype; ++i)
    {
        //SELECT MUON:
        if (type[i]==5 || type[i]==6) // Muon + or -
        {
        // OR DONT
        // Energy in GeV : mass GeV / c ^2 and momentum GeV / c
            momentum_xy = std::pow(px[i],2)+ std::pow(py[i],2);
        energyxy = TMath::Sqrt( momentum_xy);
        Et.push_back(energyxy);
        }
    }
    return Et;
}

void muonreader()
{
    //call file, can change to take filename as argument

    string fileone="/Users/Damien/Documents/rootstuff/1e3/dpmjet/DAT002206.root";
    string filetwo="/Users/Damien/Documents/rootstuff/1e3/qgsjet2/DAT002206.root";
string filethree="/Users/Damien/Documents/rootstuff/1e3/epos/DAT002206.root";

    
    //get energies
    vector<Double_t> e =    TotalE(filetwo);
    vector<Double_t> e2 = TotalE(filetwo);
    vector<Double_t> e3 = TotalE(filetwo);
    vector<Double_t> e4 = TotalE(filetwo);

    //vector<Double_t> trans = XYenergy(fileone);
    //vector<Double_t> trans2 = XYenergy(fileone);

    //count entries of e and trans
    //Int_t countTot = e.size();
    //Int_t XYTot = trans.size();
    //get tree info
    vector<Int_t> obslvl = GetBranchInt(fileone ,"ObservationLevel" );
    vector<Int_t> obslvl2 = GetBranchInt(filetwo ,"ObservationLevel" );
    vector<Int_t> obslvl3 = GetBranchInt(filethree ,"ObservationLevel" );
    vector<Int_t> obslvl4 = GetBranchInt(filethree ,"ObservationLevel" );
    
    
    //initialize histrogram he1
    TH1F *he1 = new TH1F ("he1" ," E_0 m QGSJET " ,20,0 ,200);
    he1->SetLineColor(2);

    TH1F *he2 = new TH1F ("he2" ," E_0 m DPMJET " ,20,0 ,200);
    he2->SetLineColor(8);
    
    
    TH1F *he3 = new TH1F ("he3" ," E_0 m EPOS " ,20,0 ,200);
    he3->SetLineColor(5);
    
    TH1F *he4 = new TH1F ("he4" ," E_0 m EPOS " ,20,0 ,200);
    he4->SetLineColor(7);
    
    
    //TH1I *hexy = new TH1I ("hexy" ," E_xy QGJET" ,1000 ,10 ,10000);
    //TH1I *hexy2 = new TH1I ("hexy2" ," E_xy DPMJET " ,1000 ,10 ,10000);
    //hexy->SetLineColor(4);
    //hexy2->SetLineColor(6);
    int size1 = obslvl.size();
    int size2 = obslvl2.size();
    int iterator;
    if (size1>size2) {
        iterator = size1;
    }
    else
    {
        iterator = size2;
    }
    for(long k=0; k<iterator; k++)

    {
        he1->Fill(e[k]);
        he2->Fill(e2[k]);
        he3->Fill(e3[k]);
        he4->Fill(e4[k]);
        
        //hexy->Fill(trans[k]);
    }

    
    //Chi2test doesn't work yet
    //vector<Double_t> res;
    //he1->Chi2Test(hexy,"UU",res);
    
    //make actual histo. formatting goes here.
    TCanvas *canvas_1 = new TCanvas("canvas_1", "canvas_1", 0,0,700,500);
    canvas_1->Divide(2,2,0,0);
    canvas_1 -> cd(1);
    gPad->SetLogy();
    gPad->SetLogx();
    //he1->Divide(he1,he1,1,1,option="B");
    //Double_t norm = he1->GetEntries();
    //he1->Scale(1/norm);
    
    

    //hexy->Draw("same");
    string option;
    //he1->Divide(he1,he1);
    //he2->Divide(he2,he1);
    //he3->Divide(he3,he1);
    //he4->Divide(he4,he1);
    
    //supposedly scales histo
    //Double_t norm2 = he2->GetEntries();
    //he2->Scale(1/norm2);
    he1->Draw("same");
    he2->Draw("same");
    he3->Draw("same");
    he4->Draw("same");
    
    //hexy2->Draw("same");
    
    canvas_1 -> Update();

    

}

//old tree selector
//TTree *T = (TTree*)f.Get("particle");
//T->Print();
//T->Draw("Px:Py","ParticleID == 5 || ParticleID == 6");

// SELECT MUON USING
// if (id[k]==5 || id[k]==6)
// SELECT ELECTRON USING
// if (id[k]==3)
// SELECT GAMMA USING
// if (id[k]==1)