#include "set_variables_in_root_trees.h" //you want this file, because it sets all necessary variables for you...
#include "TCanvas.h"

void simple_analysis(){
    std::cout <<"  "<< std::endl;
    std::cout <<"Load data file and set up root tree..."<< std::endl;
    
    //Get the data:
    //------------------------------------------------------
    TString file_to_analyze = "pp_data_set/pp_ppX_data.root"; //Name and directory where the root file is located... //Change this if you want to analyze a different file in a different directory
    
    //Now open the root file:
    TFile *myFile = TFile::Open(file_to_analyze);
    
    //Get the root tree:
    TTree *myTree = (TTree*)myFile->Get("data_tree"); //You need to change this IF you are analyzin a different root file
    
    //Now we are using the function defined in: 'set_variables_in_root_trees.h'
    set_tree_variables(myTree); //This function does (unfortunately) not work with root files other than 'pp_ppX_data.root'
    //------------------------------------------------------
    
    std::cout <<"...done!"<< std::endl;
    std::cout <<"  "<< endl;
    
    
    std::cout <<"Define histograms..."<< std::endl;
    
    //Define a few histograms:
    //------------------------------------------------------
    //Missing mass^2 of the beam, target, p1, p2, pi+ ,g1, g2 system...
    TH1F *mm2_2p_pip_2g = new TH1F("mm2_2p_pip_2g","",100,-0.2,0.2);
    mm2_2p_pip_2g->GetXaxis()->SetTitle("M^{2}_{x}(p_{1},p_{2},#pi^{+},#gamma_{1},#gamma_{2}) [GeV/c^{2}]");
    mm2_2p_pip_2g->GetYaxis()->SetTitle("Entries [a.u.]");
    
    //Probability of the kinematic fit including the eta-mass-constraint:
    TH1F *prob_kfit_eta = new TH1F("prob_kfit_eta","",100,0.0,1.01);
    prob_kfit_eta->GetXaxis()->SetTitle("Probability");
    prob_kfit_eta->GetYaxis()->SetTitle("Entries [a.u.]");
    
    //Looking at the kinematically fitted pi+ pi- g1 g2 invariant mass:
    TH1F *im_2pi2g_kfit_eta = new TH1F("im_2pi2g_kfit_eta","",100,0.0,1.0);
    im_2pi2g_kfit_eta->GetXaxis()->SetTitle("M(#pi^{+},#pi^{-},#gamma_{1},#gamma_{2}) [GeV/c^{2}]");
    im_2pi2g_kfit_eta->GetYaxis()->SetTitle("Entries [a.u.]");
    
    //Pi- momentum pull for the above fitter:
    TH1F *pim_mom_pull_kfit_eta = new TH1F("pim_mom_pull_kfit_eta","",100,-10.0,10.0);
    pim_mom_pull_kfit_eta->GetXaxis()->SetTitle("#pi^{-} momentum pull");
    pim_mom_pull_kfit_eta->GetYaxis()->SetTitle("Entries [a.u.]");
    //------------------------------------------------------
    
    std::cout <<"...done!"<< std::endl;
    std::cout <<"  "<< endl;
    
    std::cout <<"Loop over the tree and fill the histograms (this might take a while)..."<< std::endl;
    
    //Loop of the tree and fill the histograms:
    //------------------------------------------------------
    for(int i=0;i<myTree->GetEntries();i++){
        myTree->GetEntry(i);
        
        //Set the beam and target 4-vectors:
        set_beam_target_4vector(beam_momentum);
        
        //Calculate some masses:
        //1.) beam, target, p1, p2, pi+, g1, g2 missing mass:
        TLorentzVector mm2_2p_pip_2g_vector = beam + target - (*rec_proton1 + *rec_proton2 + *rec_pip + *rec_photon1 + *rec_photon2);
        mm2_2p_pip_2g->Fill(mm2_2p_pip_2g_vector.M2());
        
        //2.)fitted invariant mass of pi+ ,pi- , g1, g2, from the eta-mass constr. fitter:
        TLorentzVector im_2pi2g_kfit_eta_vector = (*fit_2PEta_pip + *fit_2PEta_pim + *fit_2PEta_photon1 + *fit_2PEta_photon2);
        im_2pi2g_kfit_eta->Fill(im_2pi2g_kfit_eta_vector.M());
        
        //Get the kinematic fit probability for the eta-mass constr. fitter:
        prob_kfit_eta->Fill(prob_2PEta);
        
        //Get the pi- momentum pull:
        pim_mom_pull_kfit_eta->Fill(pull_2PEta_pim[0]);
    }
    //------------------------------------------------------
    
    std::cout <<"...done!"<< std::endl;
    std::cout <<"  "<< endl;
    
    
    std::cout <<"Display the results..."<< std::endl;
    
    //Finally, display the results:
    TCanvas *c = new TCanvas("c","",1200,800);
    c->Divide(2,2);
    c->cd(1);
    mm2_2p_pip_2g->Draw();
    c->cd(2);
    prob_kfit_eta->Draw();
    gPad->SetLogy();
    c->cd(3);
    im_2pi2g_kfit_eta->Draw();
    c->cd(4);
    pim_mom_pull_kfit_eta->Draw();
    c->cd();
    
    std::cout <<"...done! Have a great day!"<< std::endl;
    std::cout <<"  "<< endl;
    
}