#include "TFile.h"
#include "TTree.h"
#include "TBrowser.h"
#include "TH2.h"
#include "TRandom.h"
#include "TGraph.h"

class SecondariesClass
{public:
	
	Double_t xPre_S, yPre_S, zPre_S, xPost_S, yPost_S, zPost_S, Ekin_Pre_S, Ekin_Post_S, Edep_S, Elost_S, ELaBr3_S, StepLength_S, xScatter_S, yScatter_S, zScatter_S , ColLength_S;
	Int_t GammaNum_S, particle_PDG_S, parent_ID_S, track_ID_S, stepNo_S, Scatter_S , user_PreVolumeIndex_S , user_PostVolumeIndex_S;
	
};

void treeSec()
{

	SecondariesClass Secondaries;

  
// file to read  
	TFile *fNU1 = new TFile("./Output1_1.root");
 

// Tree to read
	TTree *Secondaries_data = (TTree*)fNU1->Get("SecondariesTree"); 

	Secondaries_data->SetBranchAddress("SecondariesBranch",&Secondaries);

// Define histograms     
	TH1D *Energy_h = new TH1D("Energy_h","Energy Distribution",500,0,1000);
	TH3D *LaBr3Pos_h = new TH3D("LaBr3Pos_h","Gamma Distribution",500,-3,3,500,-3,3,500,5,9);   
  
	int a = -1; //variable to compare GammaNum_S to
	int b = 0;
	double ELaBr3 = 0;
 
	Int_t entries = (Int_t)Secondaries_data->GetEntries();
	cout << "Number of entries in SecondariesTree = " << entries << endl;  
	for (Int_t i = 0; i<entries; i++) { 
		 
		Secondaries_data->GetEntry(i);  
	
		// Energy deposit
		if (i%1000==0){
			 cout << "Processing entry No : " << i << endl;
			cout << "GammaNum_S = " << Secondaries.GammaNum_S << endl;
			cout << "xpre_S = " << Secondaries.xPre_S << endl;
		}
		
		if(i!=0 && i!=entries-1){
		
			if(a==Secondaries.GammaNum_S){
				ELaBr3=Secondaries.ELaBr3_S;
			}
			
			if(a!=Secondaries.GammaNum_S){
				Energy_h->Fill(ELaBr3,1.);
				a=Secondaries.GammaNum_S;
				ELaBr3 = Secondaries.ELaBr3_S;
				//cout << "zPre_S = " << Secondaries.zPre_S << endl;
				//LaBr3Pos_h->Fill(Secondaries.xPre_S,Secondaries.yPre_S,Secondaries.zPre_S);
				b++;
			}
		}
		else{
			ELaBr3=Secondaries.ELaBr3_S;
			if(i==0){
				a=Secondaries.GammaNum_S;
				//cout << "zPre_S = " << Secondaries.zPre_S << endl;
				//LaBr3Pos_h->Fill(Secondaries.xPre_S,Secondaries.yPre_S,Secondaries.zPre_S);
			}
			if(i==entries-1){
				Energy_h->Fill(ELaBr3,1.);
				b++;
			}
		}
		
		//CrystalEntrance
		
		if(Secondaries.particle_PDG_S==22 ){
			if(Secondaries.xPre_S==2.5 || Secondaries.xPre_S==-2.5 || Secondaries.yPre_S==2.5 || Secondaries.yPre_S==-2.5 || Secondaries.zPre_S==5.4){
				LaBr3Pos_h->Fill(Secondaries.xPre_S,Secondaries.yPre_S,Secondaries.zPre_S);
				cout << "IF:zPre_S = " << Secondaries.zPre_S << endl;
			}
		}
	}
	cout << "Number of primary gammas = " << b << endl;

		
  /*
// Drawing histograms 
  
   TCanvas *c3=new TCanvas("c3","Energy profile",20,10,500,500);
   c3->Divide(2,1,0.1,0.1);
   TPad* c3_1 = (TPad*)(c3->GetPrimitive("c3_1"));
   TPad* c3_2 = (TPad*)(c3->GetPrimitive("c3_2"));
   c3_1 -> Draw();
   c3_2 -> Draw();
   c3->cd(1); //first histogram
   gStyle->SetFillColor(0) ; 
   gStyle->SetOptStat(1);
   Energy_h->SetLineColor(kBlue);
   Energy_h->GetXaxis()->SetTitle("Energy");
   Energy_h->GetYaxis()->SetTitle("N");
   Energy_h->SetLineWidth(2);
   Energy_h->Draw();
   
   c3->cd(2); //second histogram
   //TCanvas *c4=new TCanvas("c4","Gamma Distribution",20,10,500,500);
   LaBr3Pos_h->SetLineColor(kBlue);
   LaBr3Pos_h->GetXaxis()->SetTitle("xPos");
   LaBr3Pos_h->GetYaxis()->SetTitle("yPos");
   LaBr3Pos_h->GetZaxis()->SetTitle("zPos");
   LaBr3Pos_h->SetLineWidth(2);
   LaBr3Pos_h->Draw("BOX");*/

   }   
   
void RootFileRead() {
   treeSec();
}