// *** Macro for drawing plots from Analysis *** //
#include <iostream>
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TCanvas.h"
#include "TMath.h"
void psi4160at8GeV10000eventsDraw()
{
		Float_t vtxprob, mom, ene, x, y, z;
		Float_t x2[4824];
		Float_t y2[4824];
		Float_t z2[4824];
		Float_t mom2[4824];
		Float_t ene2[4824];
		TBranch *b_x2;
		TBranch *b_y2;
		TBranch *b_z2;
		TBranch *b_mom2;
		TBranch *b_ene2;
		TH1D *plot1 = new TH1D("","", 200, 0, 0.2);
		TH1D *plot2 = new TH1D("","", 100, 0, 5);
		TH1D *plot3 = new TH1D("","", 200, 0, 0.2);
		TH1D *plot4 = new TH1D("","", 100, 0, 5);
		TH2D *plot5 = new TH2D("","", 1000, 0, 10, 1000, 0, 1);
		TH2D *plot6 = new TH2D("","", 1000, 0, 10, 1000, 0, 1);	
		TH1D *plot7 = new TH1D("","", 300, 0, 0.3);
		TH1D *plot8 = new TH1D("","", 300, 0, 0.3);
		TFile *thefile1 = TFile::Open("/home/apostolou/pandarootmarch14/macro/run/psi4160_2.root"); // open file with ntuples containing information
		TNtuple *thetuple1 = (TNtuple *)thefile1->Get("nd0"); // get ntuple
		TNtuple *thetuple2 = (TNtuple *)thefile1->Get("nmc"); // get ntuple
		thetuple1->SetBranchAddress("vtxprob",&vtxprob);
		thetuple1->SetBranchAddress("d0p",&mom);
		thetuple1->SetBranchAddress("d0e",&ene);
		thetuple1->SetBranchAddress("vtxvx",&x);
		thetuple1->SetBranchAddress("vtxvy",&y);
		thetuple1->SetBranchAddress("vtxvz",&z);
		thetuple2->SetBranchAddress("p", mom2, &b_mom2);
		thetuple2->SetBranchAddress("e", ene2, &b_ene2);
		thetuple2->SetBranchAddress("x", x2, &b_x2);
		thetuple2->SetBranchAddress("y", y2, &b_y2);
		thetuple2->SetBranchAddress("z", z2, &b_z2);

				
		for(Int_t j =0; j<thetuple1->GetEntries();j++)
		{
						
			thetuple1->GetEntry(j);
			Float_t s = (TMath::Sqrt(x*x+y*y+z*z));
			Float_t beta = (mom/ene);
			Float_t gamma = (1/(TMath::Sqrt(1-(beta*beta))));
			Float_t t = (s/(beta*gamma));
			Float_t bg = (beta*gamma);
			if(vtxprob > 0.01)
			{			
			plot1->Fill(t);
			plot2->Fill(bg);
			plot5->Fill(ene,beta);
			plot7->Fill(s);
			}
			
		}		
		for(Int_t j =0; j<thetuple2->GetEntries();j++)
		{
			thetuple2->GetEntry(j);
			Float_t s2= (TMath::Sqrt(x2[3]*x2[3]+y2[3]*y2[3]+z2[3]*z2[3]));
			Float_t beta2 = (mom2[1]/ene2[1]);
			Float_t gamma2 = (1/(TMath::Sqrt(1-(beta2*beta2))));
			Float_t t2 = (s2/(beta2*gamma2));
			Float_t bg2 = (beta2*gamma2);
			plot3->Fill(t2);
			plot4->Fill(bg2);
			plot6->Fill(ene2[1],beta2);
			plot8->Fill(s2);
			
		}		
	
		TCanvas *canv1 = new TCanvas("canv1","canv1");
		canv1->Divide(1,2);
		canv1->cd(1);
		plot1->Fit("expo","","",0.01,0.06);		
		plot1->Draw();
		canv1->cd(2);
		plot2->Draw();
		TCanvas *canv2 = new TCanvas("canv2","canv2");
		canv2->Divide(1,2);
		canv2->cd(1);
		plot3->Fit("expo","","",0.01,0.06);		
		plot3->Draw();
		canv2->cd(2);
		plot4->Draw();
		TCanvas *canv3 = new TCanvas("canv3","canv3");
		canv3->Divide(1,2);
		canv3->cd(1);
		plot5->Draw("BOX");
		canv3->cd(2);
		plot6->Draw("BOX");
		TCanvas *canv4 = new TCanvas("canv4","canv4");
		canv4->Divide(1,2);
		canv4->cd(1);
		plot7->Draw();
		canv4->cd(2);
		plot8->Draw();
		
}